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ABSTRACT 

The computation of theoretical pulsar populations has been a major component of 
pulsar studies since the 1970s. However, the majority of pulsar population synthesis 
has only regarded isolated pulsar evolution. Those that have examined pulsar evolu- 
tion within binary systems tend to either treat binary evolution poorly or evolve the 
pulsar population in an ad-hoc manner. Thus no complete and direct comparison with 
observations of the pulsar population within the Galactic disk has been possible to 
date. Described here is the first component of what will be a complete synthetic pulsar 
population survey code. This component is used to evolve both isolated and binary 
pulsars. Synthetic observational surveys can then be performed on this population 
for a variety of radio telescopes. The final tool used for completing this work will be 
a code comprised of three components: stellar/binary evolution, Galactic kinematics 
and survey selection effects. Results provided here support the need for further (ap- 
parent) pulsar magnetic field decay during accretion, while they conversely suggest 
the need for a re-evaluation of the assumed typical MSP formation process. Results 
also focus on reproducing the observed PP diagram for Galactic pulsars and how this 
precludes short timescales for standard pulsar exponential magnetic field decay. Fi- 
nally, comparisons of bulk pulsar population characteristics are made to observations 
displaying the predictive power of this code, while we also show that under standard 
binary evolutionary assumption binary pulsars may accrete much mass. 

Key words: binaries: close - stars: evolution ~ stars: pulsar - stars: neutron - Galaxy: 
stellar content 



1 INTRODUCTION 

Since the tentative suggestion that neutron stars (NSs) form 
from violent supernova (SN) events (Baade & Zwicky 1934a, 
1934b) and the discovery of pulsars by Hewish et al. (1968) 
the number of observed pulsars has risen dramatically. Sur- 
veys for radio pulsars have discovered over 1500 objects 
including a rich harvest of binary and millisecond pulsars 
(Manchester et al. 2005; Burgay et al. 2006). Precision tim- 
ing of pulsars in binary systems has not only allowed precise 
tests of theories of relativistic gravity (e.g. van Straten et 
al. 2001), but also given insights into their masses and the 
nature of the binary systems they inhabit (for example Ver- 
biest et al. 2007; Bell, Bailes & Bessell, 1993). Of the first 
300 or so pulsars to be discovered, only 3 were members 
of binary systems, despite the progenitor population having 
an extremely high binary fraction (> 50%; Duquennoy & 
Mayor, 1991). This paucity of binary pulsars was an impor- 
tant clue about the origin and evolution of pulsars. Clearly 
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something about pulsar generation was contributing to the 
disruption of binary systems. We now believe that the super- 
novae in which pulsars are produced impart significant kicks 
to the pulsars, that make their survival prospects within bi- 
nary systems bleak. 

Today, there are well over 100 binary pulsars known, 
and their spin periods and inferred magnetic field strengths 
offer the opportunity to attempt models of binary evolution 
and pulsar spin-up that explain their distribution in the pul- 
sar magnetic field-spin period (Bs-P) diagram. To do this 
properly, one should take models of an initial population of 
zero-age main-sequence (ZAMS) single stars and binaries, 
trace the binary and stellar evolution, including neutron star 
spin-up effects, calculate their Galactic trajectories and ini- 
tial distribution in the Galaxy, and then perform synthetic 
surveys assuming a pulsar luminosity and beaming function. 
This is what we wish to achieve. The large number of as- 
sumptions that we require to complete this effort caution 
against the absolute predictive power of such a model. For 
instance, it is easy to demonstrate that trying to use such 
a model to predict something like the merger rate of NS- 
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black hole (BH), or double pulsars/NSs based solely on a 
consideration of the total number and mass distributions of 
un-evolved ZAMS binaries would be folly. However, the rel- 
ative numbers of two or more populations can often only 
depend upon relatively few model assumptions (as shown in 
HTP02; Belczynski, Kalogera & Bulik 2002; O'Shaughnessy 
et al. 2008) . With enough observables in time we might hope 
to build up a self-consistent theory of binary and pulsar evo- 
lution. This paper, the first in a series, aims to address the 
above suggested binary and stellar pulsar evolution popu- 
lation synthesis component. Later work will combine this 
product with the kinematic and selection effect components, 
facilitating direct comparison of theory with observations. 
Therefore this paper is only a first step towards such a com- 
plete description but one that, as we will show, can already 
constrain models of neutron star magnetic field evolution 
and spin-up. 

As alluded to above, observations of pulsars take the 
form of spin measurements - both the spin period P and 
spin period derivative P - in which the surface magnetic 
field Bs and characteristic age of the pulsars are inferred 
(details are discussed when we introduce our model of pul- 
sar evolution). It is instructive to plot P vs P (hereafter 
PP) and -Bs vs P diagrams and any theoretical model must 
be able to reproduce these if it is to be successful. In Fig- 
ure [1] we show the Bs vs P diagram of ~ 1400 pulsars taken 
from the ATNF Pulsar CatalogusQ (Manchester et al. 2005). 
We see three distinct regions of the parameter space being 
populated. A large island of relatively slow spinning pul- 
sars with high surface magnetic fields (and thus high Ps) is 
joined via a thin bridge of pulsars to another, smaller, island 
of relatively rapid rotators with low magnetic fields. 

The present theoretical explanation for the distinct 
groups of observed pulsars is rotating magnetic neutron stars 
that are either isolated or reside in a binary system (see 
Wheeler 1966; Pacini 1967; Gold, 1968; Ostriker & Gunn 
1969; Gunn & Ostriker 1970; Goldreich & Julian 1969, van 
den Heuvel 1984, Colpi et al. 2001; Bhattacharya 2002; 
Harding & Lai 2006). It is those NSs that evolve within 
a binary system that may attain the low surface magnetic 
fields and low rotational periods, i.e. rapid rotators, as men- 
tioned above. The double pulsar binary J0737-3039A & B 
(Burgay et al. 2003; Lorimer 2004; Lyne et al. 2004; Dewi & 
van den Heuvel 2004) shows the distinction between rapid 
rotators and slow rotators clearly. Here we have a binary 
system consisting of a millisecond pulsar (MSP: P — 22.7 
ms and Bb = 7 x l(f G) and its 'standard' pulsar com- 
panion (P = 2.77 s and Bs = 6 X 10^^ G). Ahhough it 
seems likely that it is the process of accretion onto the NS 
that induces NS magnetic field decay (or apparent mag- 
netic field decay) this evolutionary phase is still highly con- 
tentious and theories abound on how the decay may occur. 
One such theoretical argument is of accretion-induced field 
decay via ohmic dissipation of the accreting NSs crustal cur- 
rents. This is due to the heating of the crust which in turn 
increases the resistance in the crust (Konar & Bhattacharya 
1997, 1999a, 1999b; Geppert & Urpin 1994; Urpin & Gep- 
pert 1995; Romani 1990; Urpin, Geppert & Konenkov 1997). 
An alternative explanation of accretion-induced field decay 
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Figure 1. Magnetic field vs spin period of observed pulsars 
within the Galaxy (stars, including some within globular clus- 
ters). The observations are overlaid with a cartoon depiction 
of the reason for the particulars of the pulsar parameter space 
distribution. The area covered by vertical bars primarily arises 
from the canonical pulsar birth properties and spin-down due 
to magneto-rotational energy losses. The horizontal barred re- 
gions are formed for the most part by deviation of pulsar birth 
properties from the average values from whence pulsar spin-down 
evolution commences. Forward leaning horizontal bars depict the 
region in which a luminosity law or loss of obliquity of beam direc- 
tion (or a combination of these) induces a decrease of numbers in 
the observed population. The final observed pulsar region arises 
from binary evolution, although a number of these systems are 
isolated it is possible they were formed in binaries which disrupted 
due to the explosion or ablation of the secondary star. 

consists of screening or burying the magnetic field with the 
accreted material (Lovelace, Romanova & Bisnovatyi-Kogan 
2005; Konar & Choudhuri 2004; Choudhuri & Konar 2002; 
Gumming, Zweibel & Bildsten 2001; Melatos & Phinney 
2001; Bisnovatyi-Kogan & Komberg 1974; Taam & van den 
Heuvel 1986). While yet another argument considers vortex- 
fiuxoid (neutron-proton) interactions. Here the neutron vor- 
tices latch on and then drag the proton vortices (which bear 
the magnetic field) radially, the radial direction is induced 
by either spin-up (inwards) or spin-down (outwards) of the 
NS (Jahan-Miri 2000; Muslimov & Tsygan 1985; Ruderman 
1991a, b & c). 

Below is a general overview of pulsar evolutionary 
paths. A greater level of detail is given within the very in- 
formative reviews presented by Colpi et al. (2001), Bhat- 
tacharya (2002), Choudhuri & Konar (2004), Payne & 
Melatos (2004), Harding & Lai (2006) and references 
therein. There are a number of possible binary and stellar 
evolutionary paths one can conceive that allow the forma- 
tion of a pulsar. If the end product is an isolated 'standard' 
pulsar, the star may have always been isolated and evolved 
from a massive enough progenitor star (mass greater than 
~ 10 Mq) with no evolutionary perturbations from outside 
infiuences. However, it may have been that the progenitor 
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pulsar was bound in an orbit and depending on the orbital 
parameters when the NS formed, mass loss and the veloc- 
ity kick during the asymmetric SN event (Shklovskii 1970) 
could contrive to disrupt the binary, allowing us to observe a 
solitary pulsar. Alternatively, if the secondary star was also 
massive (but originally the less massive of the two stars) 
disruption may occur in a second SN event. The creation of 
a millisecond (re-cycled; Alpar et al. 1982) pulsar requires 
the accretion of material onto the surface of the NS at some 
point during its lifetime. For this to occur the primary star is 
first required to form a NS and the binary must survive the 
associated SN event. Then, some time later, the secondary 
star evolves to overflow its Roche-lobe and initiate a steady 
mass transfer phase. Any accretion of the donated material 
onto the NS will increase its spin angular momentum, thus 
spinning up the pulsar, potentially to a millisecond spin pe- 
riod. The companion of the resultant MSP will typically be 
a low-mass main sequence star, a white dwarf (WD: a MSP- 
WD binary) or another pulsar. Formation of a double pulsar 
binary is problematic because if more than half of the sys- 
tems mass is lost in the second SN event the binary will be 
disrupted. This is unless a kick is directed toward the vicin- 
ity of the companion MSP with just the right magnitude 
to overcome the energy change owing to mass loss but not 
too strong as to disrupt the binary. If instead the binary is 
disrupted by the SN then there would be both a solitary 
'standard' pulsar and a single MSP. One other suggested 
method for producing a single millisecond pulsar is for the 
donor star to be evaporated or ablated by the extremely ac- 
tive pulsar radiation (which may be modelled in the form of 
a wind). The pulsar is then said to be a black widow pulsar 
(van Paradijs et al. 1988). 

As mentioned above, the theories of pulsar evolution can 
be tested in a statistical manner by comparing observations 
to population synthesis results. This theoretical approach 
has been adopted previously for pulsars and other stellar 
systems (some examples are: Dewey & Cordes 1987; Bailes 
1989; Rathnasree 1993; Possenti et al. 1998; Portegies Zwart 
& Yungelson 1998; Possenti et al. 1999; Willems & Kolb 
2002; O'Shaughnessy et al. 2005; Kiel & Hurley 2006; Dai, 
Liu & Li 2006; Story, Gonthier & Harding 2007 and numer- 
ous works based on StarTrack, presented in Belczynski et 
al. 2008). The level of detail in the population synthesis cal- 
culations varies, along with the methodologies the authors 
implemented. For example, Bailes (1989) considered pulsar 
selection effects in some detail, however, only roughly con- 
sidered binary evolutionary phases and the method in which 
this affects pulsar evolution. Willems & Kolb (2002) consid- 
ered NS populations and how these affect the resultant NS 
populations, however, they were not able to directly com- 
pare with observations as no selection effects were modelled. 
In a slightly different approach (empirically based) Kim, 
Kalogera & Lorimer (2003) estimated the merger rate (via 
gravitational waves) of double NS (DNS) binaries within 
the Galaxy by selecting the physical observable DNS pulsar 
properties from appropriate distribution functions for many 
pulsars and weighting these against the observed Galactic 
disk DNS pulsars PSR B1913-H16 and PSR B1534-H2. Tak- 
ing into account selection effects of the pulsar population for 
large scale surveys and producing many pulsar models they 
were able to give confidence levels for their merger rate esti- 
mates. In comparison, we will select the initial evolutionary 
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parameters from distribution functions and evolve all stars 
forward in time from stellar birth on the ZAMS through un- 
til the current time. This gives us the ability to constrain 
many stellar and binary evolutionary features while lending 
us the flexibility to, for example, compare different popula- 
tions of stars with each other and observations. In this way 
we can aim to further constrain uncertain parameter values. 
However, even at this early stage we must place a strong 
word of caution regarding the issue of parameter variation 
- strong degeneracies can apply amongst parameters (e.g. 
O'Shaughnessy et al. 2005), where a succession of parameter 
changes can mask the effect of another, so extensive mod- 
elling is required before definite conclusions can be made. 

Our goal is to create a generic code for producing syn- 
thetic Galactic populations. This will comprise three mod- 
ules: BINPOP, BINKIN and BINSFX (see the flow chart in Fig- 
ure [2] for a representation of how these fit together). The 
first module, binpop, covers the stellar and binary evolu- 
tion aspects and as such is a traditional population synthesis 
code in its own right. The second module, binkin, follows 
the positions of both binary systems and single stars within 
the Galactic gravitational potential. The third module im- 
poses selection effects on the simulations, thus giving simu- 
lated data that can compare directly to observations. This 
last consideration is in some regards the most important as 
without detailed modelling of selection effects any compar- 
ison of population synthesis simulations to observations is 
crude (Kalogera & Webbink 1998). This paper focuses on 
a description of the binpop module and, in particular, the 
pulsar population that it produces. We consider in depth 
modelling of pulsar evolution in terms of the spin period 
and the magnetic field coupled with mass accretion. Sim- 
ulating these processes will help in constraining the stellar 
birth properties of NSs and also lead to a greater under- 
standing of aspects of NS evolution such as the formation 
event itself - the supernova. It will also allow predictions of 
the composition of the Galactic pulsar population. Follow- 
up work will discus binkin with a focus on the kinematic 
evolution of the pulsar distribution within the Galaxy, and 
SN veohcty kicks - and binsfx. 

This paper is organised as follows. Section 2 gives an 
overview of the rapid binary evolution population synthesis 
code used in this research. This is followed in Section 3 by a 
detailed description of the pulsar modelling techniques that 
have been added to this code to create binpop. Section 4 
gives examples of the pulsar evolutionary pathways that can 
be followed with binpop and how these can be affected by 
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choices in the algorithm. Population synthesis results in the 
form of PP comparisons are given in Section 5 followed by 
a discussion in Section 6. In particular we wish to draw the 
readers attention to the results shown in Section FS.S.SI which 
extend beyond the basic PP description. 



2 RAPID BINARY EVOLUTION & 
POPULATION SYNTHESIS 

The research presented here makes use of the first mod- 
ule of our synthetic Galactic population code. This is called 
BINPOP and wraps the Hurley, Tout & Pols (2002: HTP02) 
BINARY STELLAR EVOLUTION (bse) cod^ (with updates as 
described in Section 3) within a population synthesis pack- 
age. 

The BSE algorithm is described in detail by HTP02 and 
an overview is given in Kiel & Hurley (2006). The aim of 
BSE is to allow rapid and robust, yet relatively detailed, 
binary evolution based on the most up-to-date prescrip- 
tions/theories for the various physical processes and scenar- 
ios that are involved. In its most basic form the BSE algo- 
rithm can be thought of as evolving two stars forward in 
time - according to the single star evolution (sse) pre- 
scription described in Hurley, Pols & Tout (2000) - while up- 
dating the orbital parameters. After each time-step the algo- 
rithm checks whether either star has over-flowed its Roche- 
lobe and depending on the result the system is evolved ac- 
cordingly, i.e. as a detached, semi-detached, or contact bi- 
nary. During these phases the total angular momentum of 
the system is conserved while orbital and spin changes ow- 
ing to tides, mass/radius variations, magnetic braking and 
gravitational radiation are modelled. 

Within BSE an effort is made to model all relevant stellar 
and binary evolutionary processes, such as mass transfer and 
common-envelope (CE) evolution. Invariably this involves 
making assumptions about how best to deal with elements 
of the evolution that are uncertain. An example that is rele- 
vant to pulsar evolution is the choices made for the nature of 
the star produced as a result of the coalescence of two stars, 
where at least one NS is involved. If the merger involves a 
NS and a non-degenerate star then a Thorne-Zytkow object 
(TZO: Thorne & Zytkow 1977) is created. Detailed mod- 
elling of these objects must deal with neutrino physics and 
processes such as hypercritical accretion and currently the 
final outcome is uncertain (Fryer, Benz & Herant 1996; Pod- 
siadlowski 1996). Possibilities include rapid ejection of the 
envelope (the non-degenerate star) to leave a single NS that 
has not accreted any mass or collapse of the merged object 
to a BH. In bse the former is currently invoked - an unsta- 
ble TZO. On the other hand, the coalescence of a NS with a 
degenerate companion is assumed to produce a NS with the 
combined mass of the two stars, unless the companion is a 
BH in which case the NS is absorbed into the BH. For steady 
transfer of material onto a NS - in a wind or via Roche-lobe 
overflow (RLOF) - it is assumed in BSE that the NS can 
accrete the material up to the Eddington limit (Cameron & 
Mock 1967). However, there is some uncertainty as to what 
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extent this limit applies as there may be cases where energy 
generated in excess of the limit can be removed from the sys- 
tem (e.g. Beloborodov 1998) . For this reason the Eddington 
limit is included as an option in BSE (HTP02). Another area 
of uncertainty is what happens to a massive oxygen-neon- 
magnesium WD when it accretes enough material to reach 
the Chandrasekhar limit - does it explode as a mass-less su- 
pernova or form a NS? Currently BSE follows the models of 
Nomoto & Kondo (1991) which suggest that electron cap- 
ture on Mg^* nuclei leads to an accretion-induced collapse 
(AIC) and a NS remnant (Michel 1987). In the formation of a 
NS via the AIC method BSE assumes that no velocity kick is 
imparted onto the NS while some other population synthesis 
works have assumed a small but non-zero velocity kick for 
these NSs (Ivanova et al. 2007; Ferrario & Wickramasinghe 

2007) . These examples illustrate some of the decision making 
involved with creating a prescription-based evolution algo- 
rithm and the interested reader is directed to HTP02 for a 
full description as well as a list of options included within 
the BSE code. 

Subsequent to HTP02 the following changes have been 
made to the BSE algorithm: 

• an option to calculate supernova remnant masses us- 
ing the prescription given in Belczynski, Kalogera & Bulik 
(2002) which accounts for the possibility of the fall-back of 
material during core-collapse supernovae (this is now the 
preferred option in bse); 

• an algorithm to compute the stellar envelope structure 
parameter. A, (required in CE calculations) from the results 
of detailed models (Pols, in preparation) - previously this 
was set to a constant whereas models show that it varies with 
mass and evolution age (Dewi & Tauris 2001; Podsiadlowski 
et al. 2003); 

• the addition of equations to detect, and account for, the 
existence of an accretion disk during mass transfer on to a 
compact object (following Ulrich & Burger 1976); and, 

• an option to reduce the strength of winds from helium 
stars. 

These changes are documented in greater detail by Kiel & 
Hurley (2006). Further additions to the BSE algorithm relat- 
ing to pulsar evolution will be described in the next section. 

There are also a number of areas where the BSE al- 
gorithm could be improved in the future. The integration 
scheme used for the differential equations in BSE (and for the 
equations introduced here) employs a simple Euler method 
(Press et al. 1992) with the time-step depending primarily 
on the nuclear evolution of the stars. It has been suggested 
that this will lead to inaccurate results for the orbital evo- 
lution, in particular when integrating the tidal equations 
(Belczynski et al. 2008) , and an implementation of a higher- 
order scheme will be a priority for the next revision of the 
BSE code. It will also be desirable to include a RLOF treat- 
ment for eccentric binaries and the addition of an option for 
tidally-enhanced mass-loss from giant stars close to RLOF 
(along the lines of Bonacic Marinovic, Glebbeek & Pols 

2008) . Another area of uncertainty is the strength of stellar 
winds from massive stars (Belczynski et al. 2008). This is a 
feature which can be varied within BSE but we do not utilise 
it or examine its effect on compact object formation in this 
work. 

Reffecting the variety of processes involved in binary 



evolution, and the uncertain nature of many of these, there 
are inevitably a substantial number of free parameters asso- 
ciated with a binary evolution code. Unless otherwise speci- 
fied we assume the standard choices for the BSE free param- 
eters as listed in Table 3 of HTP02. This includes setting 
QCE = 3 for the CE efficiency parameter. In addition we use 
the variable A for CE evolution and the remnant mass rela- 
tion of Belczynski et al. (2002) , as described above. We also 
use /inHc ~ 0.5 in the expression for the helium star wind 
strength (see Kiel & Hurley 2006 for details) and we impose 
the Eddington limit for accretion onto remnants. Where ap- 
plicable (Sections 4 and 5) we will take care to note our 
choices for any BSE parameters that we vary. 

Our aim is to produce a Galactic population of pul- 
sars. We do this using the statistical approach of population 
synthesis. The first step is to produce realistic initial popu- 
lations of single and binary stars. For each initial binary the 
primary mass, secondary mass and orbital separation are 
drawn at random from distributions based on observations 
(see Section rs. II for details). For single stars only the stellar 
mass is required for each star. We then set a random birth 
age for each binary (or star) and follow the evolution to a set 
observational time. By selecting those systems that evolve 
to become pulsars we can produce a P vs P diagram for the 
Galactic pulsar population - assuming (at least initially) 
that all generated pulsars are observable. This is repeated 
for a variety of models and comparisons made to observa- 
tions. This is the method we use here. However, it is possible 
to determine the initial parameters from a set grid and then 
convolve the results with the above mentioned initial distri- 
bution functions to determine if the binary contributes to 
the Galactic population. For more information on the grid 
method see Kiel & Hurley (2006). 



3 MODELLING A PULSAR 

We now describe our methods for evolving pulsars - both 
isolated and within a binary system. This requires a num- 
ber of additions to the BSE algorithm. We point out here 
that although most important phases of stellar and binary 
evolution are modelled we still draw the initial pulsar evolu- 
tionary parameters - spin period and magnetic field - from 
distributions. This is the method employed in all previous 
pulsar population synthesis simulations. The need for this 
arises because there is no complete theoretical method for 
modelling the SN event (Janka et al. 2006), and thus no 
strong link between the pre-SN star and post-SN star. As- 
pects such as the initial NS mass, however, are better deter- 
mined. Also, kick velocities imparted to NSs at birth are gen- 
erally selected from distributions that match observations of 
NS space velocities. We investigate the initial NS spin period 
and magnetic field distributions and extend previous studies 
by linking these to this kick velocity (see Sections 3.4 and 
5.5). We also note that although some previous pulsar popu- 
lation synthesis studies take into account selection effects to 
allow direct comparison with observations we do not do this 
here in any detail. As a first step we do consider beaming 
effects but mainly we leave this for future work. Our wish at 
this stage is to test that our evolution algorithm can pop- 
ulate the required regions of the observed parameter space. 
Considering direct number and birth rate comparisons with 
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observations will follow when the full three module code is 
complete. 



3.1 Previous efforts 

There have been many attempts at modelling NS evolu- 
tion. Some, as given in the introduction, try to generalise 
or parametrise the evolution to allow rapid computation 
for statistical considerations. Another method, however, is 
to perform detailed modelling of a single NS. Debate on 
NS evolution has concentrated mainly on the structure of 
NSs and the evolution of the magnetic field. Most statisti- 
cal studies have considered only those NSs that are isolated 
pulsars and here concern is, generally, on the magnetic field 
decay timescale or space velocity and magnetic field correla- 
tion. The semi-analytical studies of Gunn & Ostriker (1970), 
Phinney & Blandford (1981), Vivekanand & Narayan (1981) 
and Lyne, Manchester & Taylor (1985) found a magnetic 
field decay timescale of order a few million years. In compar- 
ison, studies by Taylor & Manchester (1977), van den Heuvel 
(1984) and StoUman (1987) found the decay timescale to 
be > 100 Myr. More recently population studies of pul- 
sars have tended towards a greater level of detail in the 
treatment of pulsar evolution and modelling of selection ef- 
fects. Again, differing studies find a variation of magnetic 
field decay time-scales. Those such as Bailes (1989), Rath- 
nasree (1993); Bhattacharya et al. (1992), Hartman et al. 
(1997), Lorimer et al. (1993), Lorimer, Bailes & Harrison 
(1997), Regimbau & de Freibas Pacheco (2000), Regimbau 
& de Freibas Pacheco (2001), Dewi, Podsiadlowski & Pols 
(2005) and Faucher-Giguere & Kaspi (2006) conclude that 
field decay occurs on timescales comparable to the age of 
old pulsars. While Gonthier et al. (2002) and Gonthier, Van 
Guilder & Harding (2004) suggest NSs evolve with shorter 
decay time-scales (~ few Myr). However, Faucher-Giguere 
& Kaspi (2006) suggest this later result is an artifact of 
the pulsar radio luminosity law assumed by Gonthier et. al. 
(2002) and Gonthier, Van Guilder & Harding (2004). The 
field decay time-scale is a parameter within our models that 
is allowed to vary. As we will here, Dewey & Cordes (1987), 
Bailes (1989), Rathnasree (1993), Dewi, Podsiadlowski & 
Pols (2005) and Dai, Liu & Li (2006) take into consideration 
pulsars that may evolve within binary systems and follow the 
orbital properties along with the pulsar spin and magnetic 
field. Previously either the treatment of binary evolution or 
pulsar evolution has been limited, although the work of Dai, 
Liu & Li (2006) utilises the BSE and a limited form of pulsar 
physics. No published work to date has fully considered the 
effect of magnetic field decay within both single and binary 
pulsar evolution in such detail as presented here. 



3.2 Single or wide binary pulsar evolution 

We first consider the evolution of a NS that does not accrete 
any material, either an isolated NS or one in a wide enough 
orbit to allow single star-like evolution. Although, in the 
latter case we note that aspects of binary evolution such as 
tidal interaction are followed in step with the pulsar spin 
evolution. We evolve a solitary pulsar by assuming a pulsar 
magnetic-braking model of the form. 
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(1) 



(Ostriker & Gunn 1969) with P [s/s], [G], P [s] and 
if ~ 6.3 X lO"'^* [R%sM~\ Here it is assumed that the 
beam is 90° to the pulsar equator. This is of course a 
simplification because the beam must be off-axis to pulse 
(Spitkovsky 2006). Other methods of pulsar spin evolution 
may be included in future iterations of the model. Some 
more examples are magnetar spin-down evolution (Hard- 
ing, Contopoulos & Kazanaz 1999) and spin-down of oblique 
(Contopoulos & Spitkovsky, 2006; Spitkovsky 2006) or per- 
pendicular rotators (Goldreich & Julian 1969). The mag- 
netic field is also assumed to evolve in time due to ohmic 
dipole decay which, in its simplest form, may be modelled 
as. 



Bs = -Bso exp 



TB 



(2) 



where Bso is the initial NS surface magnetic field, T is the 
age of the pulsar, tb is the magnetic field decay time-scale 
and tacc is the time the pulsar has spent accreting matter 
via RLOF. Allowing the magnetic field to decay is inspired 
by observations, where older isolated pulsars are shown to 
have lower magnetic fields than younger pulsars and is now 
a standard evolutionary model within the literature (Gumi 
& Ostriker 1969; StoUmen 1987). This apparent magnetic 
field decay is still observed even when detailed modelling of 
selection effects are taken into account (Faucher-Giguere & 
Kaspi 2006), though as mentioned above the timescale varies 
due to assumptions used within the modelling. To integrate 
these equations into the BSE we must follow the angular 
momentum of the NS using the magnetic-braking model of 
Equation [T] Differentiating the angular velocity with time 
and applying the chain rule we have the angular velocity 
time derivative. 



(3) 



Substituting Equation [T] into Equation [3] allows us to find 
the angular velocity rate of change (here it is a spin-down 
effect) in terms of the intrinsic features of the pulsar. 



'KoQ^B^ 



(I = -K-^B'i (4) 

(5) 
(6) 

where K2 ~ 2.5 x lO"^'^. It is Equation |6] which we use 
directly with the code, evolving the NS magnetic-braking 
angular momentum, Jmb with, 

Jmb = lME?tl, (7) 
5 

which is the angular momentum equation for a solid spher- 
ical rotator with radius R, mass M and angular velocity Q,. 
This is then removed from the total NS angular momentum, 
the system is updated (this includes the NS spin-down being 
transfered to the orbital angular momentum via tides) and 
we step forward in time. 

In the work presented here we are evolving the surface 
magnetic field, Bb of the NS. Strictly speaking if we were to 
compare our theoretically calculated magnetic field to ob- 
servations we should use the magnetic field strength at the 



light cylinder radius - that is the radius at which if the mag- 
netic field is considered a solid rotator it is rotating at the 
speed of light. When observations of the spin period and spin 
period derivative are made it is the magnetic field at this ra- 
dius that is being sampled. When we calculate the period 
derivative this assumed change of magnetic field is taken 
into account and as such we have two reasons for comparing 
our theoretical models to observations in the period-period 
derivative parameter space. Firstly, it is P that is actually 
being measured observationally and secondly the P calcula- 
tion naturally takes into account the light cylinder magnetic 
field (as it makes use of the pulsar magnetic moment). 



3.3 Pulsar evolution during accretion 

Allowing the decay of a NSs magnetic field during an ac- 
cretion event is believed to be the cause of the observed 
pulsar PP distribution (see Section [T]). Until now this has 
not been tested with any detailed modelling of binary, stel- 
lar and pulsar evolutionary physics. During any accretion 
phase angular momentum is transfered to the accreting ob- 
ject and it is spun-up. To follow this evolution we use the 
standard equations describing the system angular momen- 
tum as given in HTP02. While steady mass transfer occurs 
we do not allow any decrease in the NS rotational veloc- 
ity due to the magnetic-braking process described by Equa- 
tion 1^ Instead the NS angular momentum increases while 
the magnetic field decays exponentially with the amount of 
mass accreted. 



(8) 



where k is a scaling parameter that determines the rate of 
decay. Although there is no physical basis for Equation |5] if 
accretion onto a NS does cause - even at the very least - 
an apparent decay of the NS magnetic field then the process 
must be tied, in some manner, to the accreted mass (see Sec- 
tion [l]). A similar method of magnetic field evolution during 
accretion was first suggested by Shibazaki et al. (1989) and 
Romani (1990), with their equation of. 



Be 



Bso 



1 + 



AM ' 



(9) 



which is comparable to our Equation [8l Af* is a parameter 
that Shibazaki et al. (1989) fit to observations and found to 
be of order 10"* Mq. Similar to Shibazaki et al. (1989) we 
do not assume any particular model of field decay during 
the accretion phase (i.e. burial of field lines, see Section [1]), 
only that field decay following Equation [8] does occur. 

Typical accretion onto a NS has the in-falling matter in- 
teracting with the magnetosphere. Once the magnetic pres- 
sure dominates the gas pressure the matter is channelled 
along magnetic field lines to be accreted onto the magnetic 
poles (Pringle & Rees 1972). Not all mass transfer events 
within a binary system, however, result in accretion of ma- 
terial. One such example with much relevance here is the 
'propeller' (lUarionov & Sunyaev 1975) phase (Davidson 
& Ostriker 1973; lUarionov & Sunyaev 1975; Kundt 1976; 
Savonije & van den Heuvel 1977). Here we have a NS with 
a sufficiently large magnetic field and rotational velocity. 
The mass being transfered from the companion falls towards 
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the NS but instead of following the magnetic field lines to 
the surface the material is ejected from the system by the 
rapidly rotating magnetic field lines. The magnetic field acts 
as a propeller, changing rotational energy of the NS into ki- 
netic energy for the now outgoing matter. This phase has 
a noticeable effect on the NS. Due to the loss of rotational 
energy by the NS it spins down. This is therefore an im- 
portant evolutionary phase to consider during accretion and 
we have implemented a method for modelling it within BSE. 
We follow Jahan-Miri & Bhattachaxya (1994) who allow the 
difference between the Keplerian angular velocity, S7k, at 
the pulsar magnetic radius Rm and the co-rotation angular 
velocity, fi*, (Vdiff ~ (Rm) — fi*) to decide whether the 
accretion phase spins the NS up or down. The magnetic ra- 
dius is assumed to be half the Alfven radius, Rm, the radius 
at which the magnetic pressure balances the ram-pressure. 



Hwind = MIN (l, 0.01 (2 X 10" /-Bs) + O.Ol) 



(14) 



B'R" 



(2GM) M 

This gives 

Rm = 3.4 X lO'^J^B 



B^R" 



2/7 



(10) 



(11) 



The rate of change of the angular momentum for an accret- 
ing NS is then (cf. Equation 2 Jahan-Miri & Bhattacharya 
1994), 



Jacc = eVdiff-RMM, 



(12) 



where e is a parameter that allows for any uncertainties 
in the efficiency of coupling the magnetic field and mat- 
ter. During wind accretion e is considered unity, however, 
when Roche-lobe mass transfer occurs we allow e to vary 
slightly with the mass transfer rate, that is e ~ 10~^''/M 
(see Jahan-Miri & Bhattacharya 1994). Jahan-Miri & Bhat- 
tacharya (1994) conclude that this method is relatively ro- 
bust given much of the physics of accretion and capture flow 
are ignored. A similar though slightly more detailed method 
is used within the pulsar population synthesis work of Pos- 
senti et al. (1998). During the propeller phase we assume 
that the entire discarded matter leaves the binary system. 

Something that has not received a large amount of at- 
tention in binary NS population synthesis work is the as- 
sumed amount of angular momentum gained by a NS when 
accreting mass via a wind from its companion (see how- 
ever the informative works of Liu & Li 2006; Dai, Liu & 
Li 2006). Because the spin histories of NSs are important 
for our comparisons to observations, angular momentum ac- 
cretion is something that we feel should be examined. In 
particular this evolutionary phase may be important when 
modelling the apparent bridge of observed pulsars that con- 
nect the standard pulsar island to the recycled island. To 
examine this effect we multiply the angular momentum of 
the accreted mass by the efficiency parameter Hwind giving. 



2„ 2 

AJwind ~ — — windAA'fi?2^^2, 



(13) 



where Jwind is the angular momentum transferred in the 
wind, AM is the accreted mass, R2 is the companion stel- 
lar radius and Q.2 is the companion star angular velocity. 
Basing the variability of Hwind on the assumption that the 
NS magnetic field plays an integral part when NS angular 
momentum accretion occurs, we allow Hwind to vary as. 



By assuming the above form of Hwind we allow larger mag- 
netic fields to dominate the flow of angular momentum, 
while the simplicity depicts our lack of knowledge of the 
physical processes in such an evolutionary phase. An exam- 
ple of the uncertainty is shown in the work of Ruffert (1999) 
who estimates an angular momentum accretion efficiency 
between 0.01 and 0.7 (hence our lower limit for Hwind)- 



3.4 Initial pulsar parameter selection 

As mentioned earlier previous pulsar population synthesis 
work begin by selecting the initial period, magnetic field, 
mass, etc, of the pulsar from distributions, ignoring any ear- 
lier evolutionary stages the star may have passed through. If 
the distributions for each parameter are realistic this method 
is considered robust to first order. This is the method fol- 
lowed here for the initial spin period of the pulsar, Po, and 
the initial surface magnetic field, BsO- We select from flat 
distributions between Pomin and Pomax for the initial spin pe- 
riod and between PsOmin and PsOmax for the initial magnetic 
field. Typical values for these parameters are Pomin = 0.01 s, 

Pomax = 0.1s, PsOmin = 10^^ G aud P.Omax = 3 X lO" G, 

although these can vary (see Table [!}. 

In actual fact the exact initial pulsar spin period and 
magnetic field distributions are unknown. However, Spruit 
& Phinney (1998) have previously suggested that there is 
a connection between the three pulsar birth characteristics: 
imparted velocity, spin period and magnetic field strength. 
A parameter linked to the SN event and one in which we 
claim to have some form of observational contraint is Vkicii- 
The velocity kick is randomly selected from a Maxwellian 
distribution with a dispersion of Va (as in HTP02 and Kiel 
& Hurley 2006). With this in mind we trial a method in 
which the initial pulsar birth parameters Po and Pso are 
linked to Vi^icii. In this 'SN-link' method we have. 



Po = Pomin + Poav (liick/'Kr )"P , 

and 

PsO = PsOmin + PsOav ( Vkicli/ Vct )" 



(15) 



(16) 



where Poav is the average initial period and PsOav is the av- 
erage initial magnetic field. The exponents np and Uh allow 
us to parameterise the effect the SN event has on the initial 
pulsar period and magnetic field but we take both equal to 
one throughout this work. An assumption here is that a more 
energetic explosion would be able to impart a greater initial 
angular velocity onto the proto-NS than a lesser explosion, 
which may in turn introduce a more effective dynamo action 
increasing the initial magnetic field strength. 

We note here that due to the lack of a complete under- 
standing of SNe, we are at present not attempting to find the 
true pulsar initial period and magnetic field distributions hut 
trying to depict how modifying the initial pulsar birth region 
in the PP distribution affects the final pulsar PP distri- 
bution. Previous population synthesis models have differed 
on their preferred initial parameter distributions (compare 
Faucher-Giguere & Kaspi 2006 and Gonthier, Van Guilder 
& Harding 2004; Story, Gonthier & Harding 2007), however, 
any complete pulsar population synthesis must be able to re- 
produce the observed pulsars associated with SN remnants. 
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3.5 Electron capture supernova 

An important phase in NS evolution which we have touched 
on briefly already is the birth of NSs in SN events. Previ- 
ous numerical models of SNe suggested that convection may 
drive the formation of assymetries in SNe (Herant, Benz 
& Colgate 1992; Herant et al. 1994). More recently there 
have been exciting developments in this area driven by the 
successful modelling of SN explosions in multi-dimensional 
(spatially) hydrodynamic simulations (Mezzacappa et al. 
2007; Marek & Janka 2007; Fryer & Young 2007). These 
models show high levels of convection due to the stationary 
accretion shock instability (found numerically by Blondin & 
Mezzacappa 2003). This instability leads to an asymmetry 
in the explosion event and provides a mechanism for deliver- 
ing large SN velocity kicks. Observationally, proper motion 
studies of pulsars have detected large isolated pulsar space 
velocities of order 1000 km/s (e.g. Lyne & Lorimer 1994). 
Furthermore, mean pulsar space velocities appear to be in 
excess of the mean for normal field stars (e.g. Hobbs et al. 
2005). These results show the necessity of imparting veloc- 
ity kicks to NSs. As such most previous population synthesis 
works have assumed a Maxwellian SN velocity kick distri- 
bution with Va ~ 200 — 500 km/s (see above in Section |3]4]). 

Of late however, there has been a body of evidence that 
some NSs must have received small velocity kicks, of order 
50 km/s (Pfahl et al. 2002; Podsiadlowski et al. 2004). Ob- 
servationally Pfahl et al. (2002) found a new type of high- 
mass X-ray binary which exhibited low eccentricities and 
wide orbits, suggesting a low velocity kick imparted during 
the SN event. This is backed up by evidence of the relatively 
large amount of NSs within globular clusters compared to 
the estimated number if one considers they receive an av- 
erage velocity kick of Vo- > 200 km/s (Pfahl, Rappaport 
& Podsiadlowski 2002; and more recently in the theoreti- 
cal studies of Kuranov & Postnov 2006 and Ivanova et al. 
2007). The type of SN explosion that theoretically causes 
small velocity kicks, and is in vogue at present, is the elec- 
tron capture (EC) method. Basically this is the capture of 
electrons in the stellar core by the nuclei Mg, Na and 
^"iVe (Miyaji et al. 1980, see their Figure 8 and Table 1). 
This results in a depletion of electron pressure in the core, 
facilitating the core collapse to nuclear densities if the final 
core mass is within the mass range of 1.4 to 2.5 Mq (Nomoto 
1984; 1987). The bounce of material due to the halt of the 
collapse by the strong force thus drives the traditional SN 
event. Taking into account both single star evolution (Poe- 
larends et al. 2007) and binary evolution (Podsiadlowski et 
al. 2004) suggests that the initial mass range of stars that 
may explode via the EC SN mechanism resides somewhere 
within 6—12 Mq, with some dependence on metallicity 
and assumptions made in stellar modelling. Early work on 
the EC SN explosions suggested that it is a prompt event, 
in which any asymmetries would not have time to occur. 
Thus allowing the nascent NS a small SN kick velocity. This 
is in contrast to SN explosions from more massive initial 
mass stars (> 12 A^q) which produce rather slow explo- 
sions (many hundreds of ms; Marek & Janka 2007) . We note 
that the latest oxygen- neon- magnesium SN explosions (ini- 
tial stellar mass range of 8 — 12 Mq) by Kitaura, Janka & 
Hillebrandt (2006) show that the EC SN is not so prompt as 
previously thought (few hundred ms) and also energy yields 



are a factor of 10 less. However, these authors still find low 
recoil velocities for the nascent NS, again because hydrody- 
namic instabilities are unlikely to form. 

As detailed in Ivanova et al. (2007) there are a number 
of stellar and binary evolutionary pathways which may lead 
to EC SN. Although EC SN may occur from processes of sin- 
gle star evolution we concentrate here on a relatively newly 
recognised binary evolution example (Podsiadlowski et al. 
2004). A primary star with a zero-age-MS mass between 
~ 8 — 12 Mq that lives in the appropriate binary system 
may experience a mass transfer phase in which its convec- 
tive envelope is stripped from it before the second dredge-up 
phase on the asymptotic giant branch. This lack of second 
dredge-up leaves behind a more massive helium core than 
would otherwise have remained if the hydrogen-rich enve- 
lope had not been stripped off. Yet it is a less massive core 
than is left behind for stars > 12 Mq (see Podsiadlowski et 
al. 2004, Figure 1). Depending upon the rate of mass loss 
via a wind the helium core may then evolve through the 
relevant phases to explode as an EC SN (Podsiadlowski et 
al. 2004). The ability for a binary system to produce a he- 
lium star with a mass within the assumed range of 1.4 to 

2.5 AIq as given above is much greater than for a single star 
(Podsiadlowski et al. 2004). 

The possibility of EC SNe is included as an option in 
BSE for the following cases: 

• a giant star with a degenerate oxygen-neon-magnesium 
core that reaches the Chandrasekhar mass - for the stellar 
evolution models used in BSE this coresponds to an initial 
mass range of 6 — 8 Mq (for solar metalicity: HTP02); 

• a helium star with a mass between 1.4 Mq — 2.5 Mq 
(as set by Podsiadlowski et al. 2004 based on the work of 
Nomoto 1984; 1987). 

When a NS forms in these cases we use an electron capture 
kick distribution with Va = 20 km/s (see NS kick distribu- 
tion of Kiel & Hurley 2006). 

3.6 NS magnetic field lower limit 

Most detailed magnetic field decay models show a 'bottom' 
or lower limit value of the magnetic field in which no further 
decay is possible. Here, this lower limit is a model parameter, 
-Shot, and, as given in Zhang & Kojima (2006), should be of 
order 10^ - 10* G. We use 5 x 10''G throughout this work. 

3.7 Radio emission loss: The death line 

Generally, it is believed electron-positron pairs are the par- 
ticles that, when accelerated within the magnetic field of a 
NS, produce the observed radio emission. These pairs are 
assumed to form in the polar cap (PC) region of the NS 
(PC model: Arons & Scharlemann 1979; Harding & Mus- 
limov 1998) or in the outer NS magnetosphere beyond the 
polar caps (as in the outer gap model: Chen & Ruderman 
1993). Depending upon the chosen model it is possible to 
predict when or under what circumstances a NS will turn 
off its beam mechanism and thus stop being observable as a 
pulsar. These predictions form a region in the PP parameter 
space where the NS is not able to sustain pair production. 
The edge of this region is commonly referred to as the pul- 
sar death line. Until recently efforts within pulsar population 
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synthesis calculations for simulating the death-line have only 
considered isolated pulsars. Such models have been based on 
the work by Bhattacharya et al. (1992) who, following Ru- 
derman & Sutherland (1975), give the condition, 



log(p) = |log(P)-21og(/^JS.) 



El 

p2 



> 0.17 X 10 G/s^ 



(17) 



Zhang, Harding & Muslimov (2000) extend the above equa- 
tion to consider differing forms of particle scattering and 
acceleration methods but still only consider isolated pulsar 
evolution when constructing their pulsar death line models. 
These conditions, however, do not consider the physics in- 
volved in rapidly rotating NSs, such as MSPs, and can not 
account for the observed pulsar cut off at lower periods. 

Here wc must turn to the work of Harding, Muslimov & 
Zhang (2002, henceforth HMZ) who, when considering the 
PC model, make allowance for MSP spin periods and masses. 
In the long line of work completed in this field (Harding & 
Muslimov 1998, 2001, 2002; HMZ) it is shown that there 
arc a number of conditions which can be used to differ- 
entiate between the pair-production domain and the pair- 
production prohibited domain. These conditions rely heav- 
ily on two factors, (1) the type of radiation that forms a 
cascade of the aforementioned pair-produced particles and 
(2) the height at which this cascade develops - the pair for- 
mation front. For (1) HMZ calculate models bcised on three 
radiation dominant pair-production mechanisms. They as- 
sume the pair-produced particles may be created by curva- 
ture radiation (CR), which itself is released from accelerated 
particles forced to move along curved magnetic field lines, or 
via inverse Compton scattering (Hibschman & Arons 2001; 
HMZ), that is, soft photons (i.e. relatively low energy pho- 
tons) being scattered to high energies by accelerated high 
energy charged particles. If a small number of these soft 
photons have energies near the cyclotron frequency they will 
scatter with a much greater cross section and thus provide a 
larger number of high energy photons (the inverse Compton 
photons). Therefore, HMZ consider two inverse Compton 
models, resonant inverse Compton scattering and nonreso- 
nant inverse Compton scattering. For (2), if the pair forma- 
tion front occurs at an altitude less than the radius of the 
polar cap (defined by the strength of the magnetic field) the 
pulsar is said to be unsaturated. If the pair formation front 
is at an altitude greater than the polar cap radius the pul- 
sar is in the saturated regime. These regimes axe important 
because they determine in what type of electric field the 
pair-produced particles are accerlate in. Unsaturated pairs 
may accelerate in an increasing (in altituod) electric field 
(parallel to the magnetic field lines), while saturate pairs 
are accelerated in a constant electric field. 

In this work, for simplicity, we limit ourselves to onlj' 
consider the CR death line HMZ model. Our slightly modi- 
fied HMZ-CR death line equations are: 



P < 
log(P) = 



21 



log(P)-^log(/^|S.) 



3 19 
-^log{I) + -log{R)-9m 



(18) 
(19) 

(20) 
(21) 



--log(/)-hlllog(J?)-19.84, 



(22) 
(23) 



where, P < P^^ is the unsaturate regime and P^^ = 4.64 x 
W~'^Bb^^ is the CR saturated-unsaturated limit in seconds. 
The modification simply allows the use of solar units for M 
and R - B is in Gauss, P in seconds. A strength of the HMZ 
death lines is the fact that stellar evolutionary features such 
as mass and radius play a role in how the radio emission 
is formed. In fact, HMZ show how important binary mass 
transfer is in modifying the death line equations. Thus one 
parameter included within these equations is the NS moment 
of inertia, /, which we wish to calculate. To calculate / we 
are required to assume an equation of state and mass-radius 
relationship. 



2 



1 - 2.42 X 10" 



■^-2.9X10 



and 



J^ = 2.126 X W~^M~^^^. 



MR^ (24) 



(25) 



Here, I is taken from Tolman (1939) and as shown in Figure 6 
of Lattimer & Prakash (2001) this equation is a good fit 
to detailed models of NSs for differing equation of states 
down to ~ 10"'' Mq/Rq. HMZ show that the pulsar 
death line is particularly sensitive to the assumed equation 
of state, any future work on this area must take this into 
consideration. After assuming an equation of state and mass- 
radius relation we may still modify the effect of the HMZ 
death-line equations by varying /^SJJ,, the minimum pulsar 
efficiency required for pair production to occur, a parameter 
used by HMZ to fit their parameterised equations to detailed 
model results. Harding & Muslimov (2002) suggest upper 
and lower bounds on what value /^Sm should take for each 
model. For CR these values are 0.1 and 0.5 for lower and 
higher respectively. 



4 SOME EVOLUTION EXAMPLES 

To examine the newly implemented BSE pulsar evolution we 
now take a detailed look into the evolution of an isolated 
pulsar and a binary system containing a pulsar. This also 
allows us to investigate how sensitive pulsar evolution is to 
uncertain parameters in the binary and pulsar model. 



4.1 An isolated pulsar 

We begin by considering an isolated 20 Mq star initially 
residing on the main-sequence (MS) with solar metallicity, 
Z = 0.02. This star very quickly becomes a 2.3 Mq NS; 
spending only 8.8 Myr as a MS star, 0.016 Myr on the 
Hcrtzsprung gap (HG), skipping the first giant branch to 
ignite core helium burning (which lasts for 1.0 Myr) and 
then moving onto the first asymptotic giant branch, where 
after 0.02 Myr it becomes a NS. The NS, when born, is 
given a randomly selected spin period of Po = 2 x 10~^ s 
and a surface magnetic field of Bso = 5.0 x 10^^ G. This 
means the NS in this example has a period derivative of 
Po = 8.5 X 10~^^ s/s. The subsequent evolution of the NS 
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Figure 3. The evolution of pulsar particulars, P, Bs and P for differing initial parameter P and Bs values. There are two isolated 
pulsar evolutionary pathways depicted here, these are the dotted paths. Two binary pulsar evolutionary paths arc also shown (pluses 
and circles). See text for further details. 



spin period, magnetic field and period derivative are shown 
in Figures [3^), b) and c) respectively (the lower dotted 
curves in each plot). The evolution in the PP diagram is 
shown in Figure |3ji) . 

After 10 Gyr the stellar spin period, magnetic field and 
period derivative are 'observed' as, 0.76 s, 3.7 x 10* G and 
2.04 X 10~^^ s/s respectively. No pulsar has ever been ob- 
served with such parameters, so it is safe to assume that 
this NS, after 10 Gyr, would be beyond the pulsar death 
line. In fact, if we assume a death line of the form given 
by Equation [17] then the pulsar would die after a system 
time of 4555 Myr. This is the case if we assume tb = 1000 
Myr. However, if we assume instead that the magnetic field 



decays on a much faster time-scale, say, tb = 5 Myr, then 
the NS parameters evolve to a somewhat different config- 
uration. In this case the NS ends with P = 6.9 x 10~^ s, 
= 5 X lO'^ G and P = 3.6 x 10"^^ s/s. Here the NS mag- 
netic field reaches the lower limit of 5 x 10^ G in just over 
600 Myr which is about the same time that the pulsar would 
cross the death line. If we assume the star has a metallicity 
of 0.001 or 0.0001 then a BH is born directly. This is be- 
cause the boundary between NS and BH formation (which 
is governed by the asymptotic giant branch core mass which 
collapses to a 3 A4q remnant) varies slightly with metallic- 
ity. In terms of initial stellar mass this is about 21 Mq for 
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Figure 4. The evolution of pulsar particulars, P, Bb and P for differing evolutionary assumptions, tb (plus), k (circles) and propeller 
(cross). 



solar metallicity whereas for & Z = 0.0001 population it is 
19 Mq. 

4.2 A binary MSP 

We now consider the same NS/pulsar as before but place 
it within a binary system. The initially 20 Mq primary 
star (NS progenitor) has a 5.5 Mq companion in a circu- 
lar Porb = 32 day (a = 125 Rq) orbit. Initially the stars are 
assumed to reside on the ZAMS. Due to the greater mass of 
the primary star it evolves more quickly than its secondary 
companion, leaving the MS at the same time as its isolated 
counterpart of the previous example (as compared to the 
secondary, who leaves the MS after 540 Myr). The primary 



evolved to fill its Roche-lobe at a system time of 8.8 Myr. 
At this point it is in the HG phase of evolution and has lost 
~ 1 Mq in a stellar wind ~ the secondary accreted 0.285 Mq 
of this and was spun up accordingly. At the onset of RLOF 
mass transfer occured on a dynamical timescale and lead to 
formation of a common-envelope. During the common en- 
velope phase the MS secondary and the primary core spiral 
in toward each other with frictional heating expelling the 
CE. This phase is halted when the entire envelope is re- 
moved (~ 14 Mq) and the two stars are just 8 Rq apart. At 
this stage the MS secondary star is extremely close to over- 
flowing its Roche-lobe radius, however, the primary star - 
now being a naked helium star - continues to release mass in 
the form of a wind (which we assume has the form given in 
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Kiel & Hurley 2006) and due to conservation of angular mo- 
mentum begins to increase the distance between the stars. 
The secondary star accretes a small fraction of this wind 
material and grows to a mass of 5.815 Mq. Then, after just 
10 Myr the primary star goes SN ejecting a further 3 Mq 
instantly from the system to become a 1.7 Mq NS. An ec- 
centricity of 0.4 is induced into the orbit and the separation 
is now increased to 15 Rq. By the time the secondary star 
evolves enough to initiate steady mass transfer, at a sys- 
tem time of 67 Myr, the orbit has been circularised by tidal 
forces acting on the secondary stars' envelope. The accretion 
phase lasts for 2800 Myr. Within this time the secondary 
star evolves from the MS to the HG and then onto the first 
giant branch. During accretion the secondary loses 5.6 Mq, 
the NS accretes 1.1 Mq and the separation increases out 
to 32 Rq. Beyond this the system does not change greatly. 
After a short time (~ 30 Myr) the secondary evolves to be- 
come a 0.2 Mq helium WD (HeWD) - being stripped of 
its envelope the WD does not develop any deep mixing to 
form a higher metalicity content in its envelope (e.g. such 
as carbon-oxygen or oxygen-neon-magnesium WDs). At the 
end of the simulation we are left with a 2.885 Mq NS with 
a 0.2 Mq WD companion in an ~ 10 day orbit. This is 
typical of the parameters of observed pulsar-HeWD binaries 
(van Kerkwijk et al. 2004). We have assumed a maximum 
NS mass of 3 Mq, however this mass limit is not well con- 
strained. If we were to assume a maximum NS mass of only 
2 Mq the system would end as a WD orbiting a BH. 

In terms of pulsar evolution, the magnetic field ver- 
sus period parameter space covered by the NS is shown in 
figures [3] and |4] for a variety of parameter value changes. 
To begin we look at how simply changing the initial pe- 
riod and magnetic field changes the NS particulars over 
time. We direct the reader to Figure [3] which shows the 
clear differences in the binary pulsar life if it starts its life 
with Po = 2.5 X 10"^ s, BsO = 5 x 10^^ G and therefore 
Po = 9 X lO^^'' s/s (pluses) or if it starts its life with 
Po = 9.3 X 10"^ s, Pso = 1.4 X 10^^ G and Po = 2 X 10"" s/s 
(circles) . There is a marked difference between the spin evo- 
lution of the pulsars - the evolutionary pathways depend 
greatly on the strength of the pulsars magnetic field. It is 
obvious from these curves where accretion began and be- 
yond this the pulsar evolution is similar for the two stars. 
Comparison can also be made with the isolated pulsar evo- 
lution given in the previous example (lower dots), initially 
these two pulsars evolve similarly until accretion occurs. 
Another isolated pulsar is shown (higher dots), with Po = 
3.4 X 10"^ s, Pao = 4.1 X 10^^ G and Po = 4.9 x 10"" s/s. 
This isolated pulsar spins down much faster than the other 
three pulsars (see Figure |3^) and shows the effect a greater 
magnetic field has on the spin evolution of a pulsar. Once ac- 
cretion occurs the spin and magnetic field evolution for both 
binary pulsars changed dramatically. The NSs are eventu- 
ally spun-up and both possess sub-millisecond spin periods, 
while the magnetic fields decay rapidly and both reach their 
lowest limits (5 x 10*^ G) in ~ 30 years. The P and P param- 
eter space covered by the three pulsars is given in Figure [3J1. 

At this stage we have not considered allowing the NS 
to receive a velocity kick during the SN event. If we allow 
a large, randomly orientated velocity kick of ~ 190 km/s to 
occur during the SN then the system is disrupted. However, 
if a smaller velocity kick is given, say, ~ 50 km/s not only 



does the system survive but it passes through a complex set 
of mass transfer phases including two RLOF phases and a 
handful of symbiotic phases. The first RLOF mass transfer 
phase, which began at a system time of 75 Myr, caused the 
NS to be spun up to 0.02 s, accrete 0.04 Mq and lasted 
for 0.2 Myr. The secondary star lost 4.8 Mq of mass and 
the separation reduced to 5.8 Rq. During this phase the 
magnetic field decayed to its bottom value during the first 
Roche-lobe phase and the pulsar crosses any assumed death 
line and can not be observed beyond this evolutionary point. 

4.3 Binary evolution variations 

The examples of binary MSP formation above assumed evo- 
lutionary parameters of tb = 1000 Myr, k = 10000 and no 
propeller evolution. We now reconsider the example with 
Po = 2.5 X 10"^ s and BsO = 5 x 10^^ G and vary these 
assumed parameters one at a time. The variations are (see 
Figure lU : 

• Assume re = 5 Myr (plus symbols within Figure [Jjl; 

• Decreasing fc to (circles); 

• Allow propeller evolution to occur (crosses); 

As expected a lower value of re forces the magnetic field 
to decay rapidly while causing little spin-down of the pulsar 
spin period to occur. The binary evolution is not noticeably 
influenced due to this evolutionary change. 

Decreasing k affects the rate of magnetic field decay 
during accretion (see equation (Sjl as depicted clearly in Fig- 
ure UJd. This parameter change seems one we would defi- 
nitely regard with skepticism because it populates entirely 
the wrong region of the PP parameter space, in terms of 
what observation suggest. The binary evolution, again, does 
not change significantly. 

Generally it is thought that if there is going to be a pro- 
peller phase it will occur at the beginning of the accretion 
phase (Possenti et al. 1998), when the pulsar is rotating 
rapidly enough and the magnetic field is sufficiently large 
enough. This is what occurs here, however, unlike the pro- 
peller phase considered by Possenti et al. (1998), our method 
causes it to initiate multiple times during this evolutionary 
example. During the first propeller phase the NS is spun- 
down, reaching a spin period of 5 x 10^ s before the phase 
ends and accretion onto the NS begins. Although the RLOF 
phase lasts for ~ 2000 Myr - shorter than without a pro- 
peller phase - the accretion of mass onto the slowly rotating 
NS occurs extremely rapidly and unfortunately is only cal- 
culated by the BSE over a single time-step. This is the reason 
for the large ammount of magnetic field decay in one step. 
Modifying the time-stepping around this region does not af- 
fect the evolution greatly, however, it does slow the speed 
at which we may evolve many systems. The first propeller 
phase although very efficient in braking the NS spin was 
short lived, lasting only ~ 2.3 Myr. While beyond the first 
mass transfer phase the system moved once again into pro- 
peller evolution, this time lasting for ~ 550 Myr and spin- 
ning the NS down to a ~ 1000 s spin period. Once again, at 
a system time of ~ 580 Myr the NS accretes material and is 
spun in one time-step up to ~ 0.04 s spin period. The third 
propeller phase, initiated directly after the accretion phase, 
continues until the end of the RLOF mass transfer phase 
and manages to spin the NS to a period ~ 0.08 s. At the 
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end of the RLOF phase the companion, now on the first gi- 
ant branch and 40 Rq away fi'om the NS, loses ~ 0.001 Mq 
of matter the majority of which is collected by the NS. This 
spins up the pulsar to 0.07 s. We end with a 0.3 AIq HeWD, 
slightly greater than that created without the propeller in- 
cluded. While the NS is 1.78 Mq, slightly smaller than when 
we do not assume the propeller phase to occur. 

From the simple analysis above we can see how effective 
the propeller phase is at braking the pulsar spin. Although 
the propeller phase arrests the majority of spin-up, and in 
fact obstructs the formation of a MSP (that would otherwise 
have formed), it is an important evolutionary assumption to 
test on populations of NSs. We have also shown that both 
the accretion-induced field decay parameter k and the stan- 
dard magnetic field time-scale parameter tb are important 
parameters to regard when modeling the entire pulsar pop- 
ulation. We next consider populations of pulsars and how 
parameter changes affect the morphology of these popula- 
tions. 



5 RESULTS AND ANALYSIS 
5.1 Method 

A binary system can be described by three initial parame- 
ters: primary mass, M, secondary mass, m, and separation, 
a (or period, Porb). A fourth parameter, eccentricity, may 
also vary, however, for simplicity we assume initialy circular 
orbits. For the range of primary mass we choose a minimum 
of 5 Mq and a maximum of 80 Mq . To first order this covers 
the range of initial masses that will evolve to become NSs 
via standard stellar evolution. Stars below ~ 10 Mq will 
generally become WDs and stars above ~ 20 Mq will end 
up as BHs. However, binary interactions tend to blur these 
stellar mass limits. As sugested in HTP02 the production of 
stars with M > 80 Mq is rare when using any reasonable 
IMF, hence the choice of upper mass limit. The lower limit 
is one that should allow the production of all interesting NS 
systems to be formed, that is, all systems which have interac- 
tion between the NS and its companion. This includes stars 
with initial masses below 10 Mq that accrete mass and form 
a NS instead of a WD. Expanding the upper limit accounts 
for stars with M > 20 Mq that lose mass and become a NS 
instead of a BH. The secondary mass range is from 0.1 Mq 
to 40 Mq. The lower limit here is approximately the mini- 
mum mass of a star that will ignite hydrogen-burning on the 
MS. These mass ranges will allow the production of the ma- 
jority of NS binaries and enable us to investigate complete 
coverage of the PP diagram. For the initial orbital period 
we consider a range between 1 and 30000 days. The birth 
age is randomly selected between and 10 Gyr where the 
latter reflects an approximate age of the Galaxy. 

Instead of calculating birth rates and numbers of sys- 
tems of interest, our focus, here, is on comparison with the 
observed pulsar P vs P diagram and what this can teach 
us about uncertain parameters in pulsar evolution. These 
include parameters such as the magnetic field decay time 
constants in both the isolated and accreting regimes along 
with the effectiveness of modelling propeller systems to spin- 
down pulsars to observed periods. We also aim to test how 
different evolutionary parameter changes affect the relative 
populations of binary and isolated pulsars. 



In distributing the primary masses we use the power 
law IMF given by Kroupa, Tout & Gilmore (1993), 

r 0.035M-^-^ if 0.08 s; M < 0.5, 
$ (M) = <^ 0.019M-^-^ if 0.5 s; M s; 1.0, (26) 
[ O.OigM"^ '' if 1.0 sC M < oo. 

This is the probability that a star has mass M [Mq]. The dis- 
tribution of secondary star masses is not as well constrained 
observationally. We assume that the initial secondary and 
primary masses are closely related via a uniform distribu- 
tion of the mass-ratio (Portegies Zwart, Verbunt & Ergma 
1997; HTP02), 

^M = ^, (27) 

where the mass-ratio ranges from 0.1 /M to 1.0 The initial 
orbital period distribution we take to be flat in log (Porb). 
Completely sampling these distributions allows us the full 
range and distribution of initial parameters, so we are mod- 
elling the natural selection of stellar systems as suggested 
by observations. 

The results of six simulations (see Table [T| in the form 
of PP plots are given below (Sections 15.21 to I5.7|l . These 
are Models A through F. We begin by simulating a rather 
naive stellar population (Model A) to which we build pulsar 
evolutionary phases and assumptions upon in increasingly 
greater complexity until we feel the results best represent 
the required pulsar PP parameter space as compared to ob- 
servations. We note from the onset that we do not attempt 
to model the slow rotating high magnetic field pulsars in 
the upper-right region of the PP parameter space. As dis- 
cussed in Section [6] we leave these potential magnetars for 
future work. Modifying the assumed initial NS spin period 
and magnetic field (see Section 15. 6p quantifies the impor- 
tance of choices made in this area. For the most part the 
initial period and magnetic field ranges (given in Table [1} 
are selected so as to allow good coverage of the PP parame- 
ter space. At this stage we do not suggest that these ranges 
are completely realistic but we do wish to show that it is 
possible to produce a complete coverage of the pulsar PP 
parameter space with them. In Section [5]8] we take Model F 
and look at some further variations such as inclusion of elec- 
tron capture supernovae and beaming effects. This includes 
an examination of the pulsar primordial parameter space 
("Section 15. 8. 4|l and how changing evolutionary assumptions 
modifies the pulsar birth properties. 

Unless otherwise stated all models have the binary and 
stellar evolutionary parameters of: solar metalicity Z — 0.02; 
a maximum possible NS mass of 3 Mq and acE = 3. SN 
kicks are given to NSs, while to keep the required number 
of models down to a minimum we only use the curvature 
radiation death line model with /p^Sm — 0.15. Although we 
know the CR death line to be insufficient in terms of mod- 
elling recycled pulsars it is adequate for our purposes here 
(for more see Section [6]). For each model we evolve 10 pri- 
mordial systems, all of which are initially binary systems. 
Integrating over the Kroupa, Tout & Gilmore IMF with the 
mass limits set above, modeling 10^ primordial systems is 
roughly representative of 2 x 10^ systems in the Galaxy if the 
full mass range is allowed. One way to think of this is that 
we are assuming a stellar birth rate of 10~^ systems/year. 
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assuming the age of the Galaxy is 10 Gyr. In mass, we are 
only evolving ~ 0.1% of the Galaxy assuming it contains 
~ 1 X 10^^ Mq. This is something that in the future we will 
wish to increase. 

5.2 A starting model 

We start by demonstrating the results of pulsar population 
synthesis when no field decay owing to the accretion of ma- 
terial onto a NS is allowed. This is done by setting fc = in 
Equation [H] The result is Model A for which the pulsar evo- 
lutionary parameters are given in Table [1] and the PP distri- 
bution is shown in Figure [S] Clearly from Figure [S] there are 
a great number of rapidly rotating pulsars with large spin 
period derivatives. This is not observed which means that 
this model fails to decay the magnetic field quickly enough 
during accretion. Therefore, Model A seems physically un- 
realistic. One may be skeptical of this result based on one 
model, however, similar PP distributions occur even if we 
vary additional parameters such as Hwind (to lower the an- 
gular momentum transfered during accretion) and/or using 
lower values of tb. It is interesting to observe the different 
levels of pulsar densities in the spun-up region of Figure [5] 
- there is almost a wave-front effect towards MSP periods. 
The pulsars in these regions have differing companion types 
and orbital configurations indicating different mass transfer 
histories. Here the fastest rotating pulsars generally have MS 
companions and circular orbits. In these cases mass transfer 
occurs via RLOF. For the pulsars located in the spin-up re- 
gion truncated at P ~ 0.002 s companion types are mainly 
giant and white dwarfs and include many cases of eccen- 
tric orbits. This indicates mass transfer via a stellar wind 
from a giant star in a binary with a relatively long orbital 
period. Systems which have a shorter orbital period and ini- 
tiate RLOF from a MS companion can accrete a greater 
ammount of material and hence are spun-up further to left 
in the PP diagram. 

5.3 Magnetic field decay during accretion 

Model A shows that some form of accretion induced mag- 
netic field decay is required to suitably populate the recy- 
cled region of the PP diagram. A first attempt at mod- 
elling this physics is completed in Model B (see Table [1]) , 
where we fc = 10000 in Equation |8l The resultant PP pa- 
rameter space is shown in Figure |6l This simulation gives 
good coverage of the observed pulsar regions - the standard 
and recycled islands. However, there is an over-production 
of systems linking the two pulsar PP islands. Many of these 
pulsars are paired with early type stellar companions (pre- 
degenerate stars) with strong winds that may disrupt the 
radio beam (see Section [5.8.2[) . A number of these systems 
are also beyond what is considered the general 'spin-up' line. 
The spin-up line is a theoretical construct which suggests 
that an accreting pulsar will eventually reach some equilib- 
rium spin period and therefore end pulsar spin-up during the 
accretion phase. The spin-up line, however, is a very uncer- 
tain model and the majority of assumptions that go into it 
are highly variable (see Arzoumanian, Cordes & Wasserman 
1999 and references therein for details), and the line is sen- 
sitive to changes for many parameter values. We do not di- 
rectly impose the spin-up line in our models but instead aim 
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Figure 5. The PP diagram for Model A. The relevant parameter 
values are given in Table[T]and on the figure. Plus symbols repre- 
sent the observed pulsars, while grey dots represent each pulsar 
the simulation produces. 



to have its effect replicated by the physical processes that 
make up our binary/pulsar evolution algorithm. As such we 
would clearly want to limit the production of medium-low 
period high-P pulsars that appear in Model B. This is a mo- 
tivator for testing the affect of including propeller physics. 
One may also notice the line of MSPs at the bottom right of 
Figure [S] This line appears in all the PP figures shown here 
and is caused by the assumed limit of magnetic field decay 
(Section [32]). 



5.4 Propeller evolution 

Propeller evolution (Section 13. 3p is permitted in Model C, 
which is otherwise the same as Model B. The results are 
shown in Figure [7] which exhibits three main features that 
make modelling the propeller evolution beneficial. The first 
is the slight bump from the standard pulsar island in slower 
spin period pulsars (P > 5 s). This slight bump is directly 
caused by the propeller phase, and allows us to model those 
pulsars who reside next to the standard pulsar island. Al- 
though, we note that those pulsars could also be modelled if 
we assume that the magnetic field does not decay other than 
during accretion. The second feature of interest is the cut off 
of millisecond pulsars with period derivatives between lO"'^^ 
and 10~^*, as compared to Model B. Although this cut-off 
occurs too soon (in terms of pulsar period) it does provide 
a method for producing this observed limit in period. The 
third feature is lack of pulsars with period derivatives be- 
tween 1Q~" and 10~" and spin periods of around 0.01. 
Therefore, when considering PP parameter space, we find 
that including the propeller phases provides much better 
agreement with models that impose a spin-up line directly to 
restrict pulsars appearing above the observed pulsar bridge. 
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Table 1. Characteristics of the main set of models used in this work. The first column gives an identifying 
letter for each model. Within the second column we indicate the value of the magnetic field decay time-scale 
which is followed by the choice of accretion induced magnetic field decay parameter. The fourth and fifth 
columns indicate whether or not propeller evolution is considered and whether a link between SN strength 
and initial magnetic field and pulsar period is assumed. The next four columns indicate the maximum and 
minimum values for the initial distributions of pulsar spin period and magnetic field. We then modify the 
ammount of angular moment accreted by a NS when the companion is lossing mass in a wind, this is governed 
by the parameter S^i^d- The penultimate column indicates if the electron capture SN kick distribution is 
considered and the final column indicates whether beaming and accretion selection effects are considered. 



Model TB [Myr] fc Propeller SN link -Pominls]^ -Pomax [s] -Bomin[G] Somax [GJ^ S^ind EC SN Beaming 



A 


1000 





No 


No 


0.01 


0.1 


1 X 10^^ 


3 X 10^^ 


1 


No 


No 


B 


1000 


10000 


No 


No 


0.01 


0.1 


1 X 10^^ 


3 X 10^3 


1 


No 


No 


G 


1000 


10000 


Yes 


No 


0.01 


0.1 


1 X 10^^ 


3 X 10^^ 


1 


No 


No 


D 


5 


10000 


Yes 


No 


0.01 


0.1 


1 X 10^^ 


3 X 10^^ 


1 


No 


No 


E 


2000 


3000 


Yes 


Yos 


0.02 


0.16 


5 X 10^^ 


4 X 10^^ 


1 


No 


No 


F 


2000 


3000 


Yes 


Yos 


0.02 


0.16 


5 X 10^^ 


4 X 10^^ 


variable 


No 


No 


Fb 


2000 


3000 


Yes 


Yos 


0.02 


0.16 


5 X 10^^ 


4 X 10^^ 


variable 


Yes 


No 


Fc 


2000 


3000 


Yes 


Yes 


0.02 


0.16 


5 X 10^1 


4 X 10^^ 


variable 


Yes 


Yes 


Fd 


2000 


3000 


No 


Yos 


0.02 


0.16 


5 X 10^1 


4 X 10^^ 


variable 


No 


No 



For the SN link models this value is the average initial spin period Poav, not the minimum. 
^ For the SN link models this value is the average initial magnetic field Boavi not the maximum. 
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Figure 6. The PP diagram for Model B. The relevant parameter 
values are the same as Model A except for the accretion induced 
field decay parameter which is now fc = 10000. Plus symbols 
represent the observed pulsars, while grey dots represent each 
pulsar the simulation produces. 
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Figure 7. The PP parameter space for Model C. Plus symbols 
represent the observed pulsars, while grey dots represent each 
pulsar the simulation produces. 



5.5 Magnetic field decay timescale 

Model D differs from tiie evolutionary assumptions in Model 
C by assuming tb ~ 5 Myr (reduced from 1000 Myr). This 
rapid field decay has the effect of shifting the standard pulsar 
island to the left in the PP diagram (see Figure [8} . Even 
so, looking at Figure [8] we see it is not possible to properly 
simulate the standard pulsar island with tb = 5 Myr, unless 
a completely contrived initial period range is used. Also in 
Figure |8] it is noticable that even before consideration of 



selection effects there is a large underproduction of recycled 
pulsars. The rapid field decay causes the majority of spun- 
up pulsars to sit on the lower magnetic field limit, a region 
of the PP parameter space in which any complete death 
line model should remove. As such, assuming such a low 
value of magnetic field decay, as has been used in the past 
by Gonthier et al. (2002), for example, is unrealistic. We 
find that tb of 100 Myr or greater is required to produce a 
reasonable PP distribution. 
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Figure 8. The pulsar PP parameter space for observed and 
those simulated pulsar produced by Model D. This facilitates 
comparison between a large and small tb when compared to Fig- 
ure [7] Plus symbols represent the observed pulsars, while grey 
dots represent each pulsar the simulation produces. 

5.6 Updated initial pulsar distribution: SN link 

With the propeUer phase producing what seems to be a more 
reahstic population in Model C than that which occurred in 
Model B, while also constraining tb 5S 100 Myr in Model 
D, we now conduct a trial to see what effect linking the SN 
velocity kick magnitude with the initial period and mag- 
netic field (see Section 13. 4p has on the pulsar distribution 
in PP. To do this we make use of equations 1151 and 1161 as- 
suming rip and Wb are unity. The corresponding values of 
Poav and -Boav are given in Table [1] (see Model E). Note 
that the distribution of initial pulsar period and magnetic 
field is modified (as compared to Models A through C) to 
allow better coverage of the standard pulsar island in the 
observed PP space. The value of the lower initial magnetic 
field limit, BaOmin, was used due to suggestions that there 
are NSs which are born weakly magnetised (~ 5 x 10^^ G; 
Gotthelf & Halpern 2007) . We also allow a greater coverage 
of the slower rotating pulsars by assuming tb = 2000 Myr. 
The previous value of tb = 1000 Myr under populated the 
slower pulsars when used in conjunction with this model. 
Conversly if the value of tb is greater than 2000 Myr the 
lower region of the standard pulsar island is increasingly 
truncated - the magnetic field decays too slowly to suffi- 
ciently decrease the period derivative. These changes result 
in Model E and its PP diagram is shown in Figure E] There 
are now striking similarities between the observed and sim- 
ulated PP parameter spaces, particularly at the standard 
pulsar island. Three main features in Figure [9] of worthy 
note are as follows: (1) the production of a wall of young 
pulsars, which indicates that those pulsars born with strong 
magnetic fields also initially spin at a faster rate to those 
born with lesser magnetic fields (which is what we expect 
from Equations 1151 and I16|) : (2) the only method for pul- 



sars to have periods of ~ 4 seconds or slower, in this model, 
is to pass through a propeller phase; and (3) the observed 
group of high P (~ 10"" - 10"^^ [s/s]) but moderately 
rapid rotating pulsars (~ 0.01 — 0.2 [s]) are reasonably well 
modelled. It should be pointed out here that within Model 
E there is a different reason why the PP area described 
in point three is populated as compared to the explanation 
given for Figure [T] Within Model E these systems arise due 
to mass transfer, while the assumed formation mechanism 
of these pulsars (for Figure [ij was that these where newly 
formed NSs and born there. 

In terms of the shape of the PP parameter space, we 
are able to simulate some regions of the observations rather 
accurately. To begin we only slightly over-estimate the cut- 
off of pulsars at small spin periods. This cut-off arises from 
a combination of the assumption of k and the propeller evo- 
lutionary phase. The death line does a satisfactory job at 
limiting the semi-recycled pulsars however it fails in limit- 
ing both the slower rotating pulsars and the fully recycled 
pulsars. This failure is not unexpected because the curva- 
ture radiation death line is not a completely realistic model 
when considering the full pulsar population (see Section [S] 
for futher discussion on this point). 

What is not noticable in Figure is the number of pul- 
sars that evolve through a Thorne-Zytkow object phase. As 
explained earlier a TZO is a remnant core (NS or BH) sur- 
rounded by an envelope formed from the non-compact star 
with which it has merged. In BSE these objects are treated 
as unstable so that the envelope is ejected instantaneously. 
However, we are left with the question of what to do with 
the pulsar parameters. For now we keep this evolutionary 
phase consistent with our CE phase were we assume that 
the spin period and magnetic field are re-set by this pro- 
cess. This assumption is prefered at this stage because if we 
assume the pulsar is left unchanged we find many isolated 
pulsars rotating near their assumed break-up spin velocities 
(P ~ 4 X 10"** s) with a large range of spin period derivatives 
(P < 10^^^). We note here that our assumption regarding 
these objects is purposefully simplistic to reflect our lack 
of knowledge about what occurs during the formation and 
evolution of a TZO but other possibilities are discussed in 
Sections [5X2I and [6l 

Once again, there is an over production of pulsars above 
the observed pulsars which link the two pulsar islands. In the 
next section we again increase the detail of our model in an 
attempt to constrain this PP region. 



5.7 Testing wind angular momentum accretion, 
the H„ind parameter 

The effect of modifying Hwind (Model F) is shown in Fig- 
ure [To] and is dramatic when compared with Figure |9l This 
modification truncates the width of the pulsar bridge which 
was the aim of this model. The pulsar bridge that stretches 
between the two pulsar islands compares rather accurately 
to the observations (in terms of PP area covered). This is 
our prefered model and we now perform an in depth analy- 
sis of this model, including further parameter changes (see 
Table [1|. 
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Figure 9. The pulsar PP parameter space for Model E. Plus 
symbols represent the observed pulsars, while grey dots represent 
each pulsar the simulation produces. 
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Figure 11. The binary (circles) and isolated (stars) pulsar PP 
parameter space for Model F. 
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Figure 10. The pulsar PP parameter space for Model F. Plus 
symbols represent the observed pulsars, while grey dots represent 
each pulsar the simulation produces. 



5.8 Model F: analysis and further additions 

Before we embark on our further examination of Model F 
we must remind the reader that - at this stage - our pa- 
rameter space analysis is incomplete. Therefore, although 
we provide evidence for the formation of such systems as: 
isolated MSPs, eccentric double neutron stars and possi- 
ble existence of Galactic disk eccentric MSP-BH binaries; 



we must stress that it is possible for many realistic model 
parameter changes to modify the resultant formation prob- 
ability of these systems. However, we provide this section 
as an example of typical binary, stellar and pulsar features 
that our future complete pulsar population synthesis pack- 
age will be able to constrain when one is able to compare 
theoretical results directly with observations. This section 
also provides evidence of how some evolutionary assump- 
tions further modify the final pulsar population. 

Figure[TT]depicts the isolated and binary pulsar systems 
in the PP parameter space for Model F. Focusing on MSPs, 
the relative number of binary to isolated MSPs is 0.49. Note 
that we are restricting our comparison to those pulsars that 
have a magnetic field greater than 6 x lO'^ G and less than 
1 X 10® G (model pulsar period region is P < 0.02; see Fig- 
ure [T] for observational comparison) . In doing this we are 
ignoring those pulsars that sit on the bottom magnetic field 
limit. To examine the MSP population of Model F in fur- 
ther detail we show the initial primary and secondary mass 
distributions of the isolated and binary MSPs in Figure [T2b 
and see that these differ. Notably the average initial mass 
of the secondary stars is typically greater for binaries that 
go on to produce isolated MSPs. What this is showing is 
that the isolated MSPs are passing through two SN events 
and it is the second SN that disrupts the system. This pro- 
vides us with a simple method for producing isolated MSPs 
rather than requiring the addition of some ad-hoc model 
where the MSP is destroying the companion. This suggested 
formation mechanism of isolated MSPs is not restricted to 
Model F alone but to all models completed so far, in fact, 
due to our restriction of the accreted wind angular momen- 
tum we are now producing less isolated MSP systems than 
previously. Therefore, we believe this is an important evolu- 
tionary phase that should be placed under further scrutiny 
within the community. 
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Figure 12. The binary (circles) and isolated (dots) MSP initial 
mass parameter space for Model F in a and Model Fb in b. All 
systems here fulfill the 6 X lO'^ ^ Bs ^ 1 x 10^ G and P < 0.02 
criterion. 

A number of these disrupted MSPs arise from systems 
which contain initial secondary masses lower than the as- 
sumed upper mass limit for the EC SN model as shown in 
Figure [12] This suggests that if we aUow EC SNe and adopt 
the associated lesser Va value (see Section I3.5|l that many 
of the now isolated MSPs may retain their binary compan- 
ion. This leads us to consider the effect the combined EC 
SN and standard SN velocity kick distribution has on our 
model pulsar population. 

5.8.1 A new supernovae velocity kick distribution: 
electron capture supernovae 

We have repeated Model F with the electron capture SN 
kick distribution taken into consideration to produce Model 
Fb. We find that Models F and Fb do not differ greatly 
from each other in terms of the regions of the PP param- 
eter space covered (hence we do not show Model Fb). Sur- 
prisingly the binary to isolated MSP ratio has only risen to 
0.53 in comparison to Model F. However, Model Fb does 
contain a greater number of binary MSPs. The initial mass 
parameter range of MSPs that arise in Model Fb is shown 
in Figure 112b ) . It is obvious that many more binary MSP 
systems are now being born in and around the ~ 10 A/© 
primary mass and ~ 7 — 10 Mq secondary mass range as 
expected. Not so expected was the number of newly formed 
isolated MSPs that are also derived from this primary and 
secondary mass region. It appears that many primary mass 
stars are forming NSs via the electron capture method, al- 
lowing them to remain within a binary system and receive 
angular momentum from their companion which then also 
goes SN, disrupting the system. In particular as Figure IT^) 
shows, there are a greater number of isolated (and binary) 
MSPs being born from the 20 Mq primary mass region than 



in Model F (Figure [Hj,) . A thorough analysis of the ECSN 
model can not be completed at this stage as considerations 
of NS space velocities are required for any complete exami- 
nation of this assumption. However, we have demonstrated 
some effects the ECSN model imparts on the pulsar popu- 
lation. As such it seems likely that the EC assumption will 
play an important role in future NS population synthesis 
works. Especially if comparison to observations of orbital 
period, eccentricity and space velocities are attempted. 



5.8.2 Beaming effect and accretion: a start on selection 
effects 

The source of radio emission from a pulsar occurs at the 
magnetic poles, otherwise known as the polar caps. The ra- 
dio waves are well collimated and form a beam. This beam 
sweeps across the sky (as observed from the pulsar) and may 
be observed if the beam crosses the observer line of sight. 
Some regions of the pulsar sky are not covered by the beam 
- the beaming effect considers the chance that a particular 
pulsar is beaming across our line of sight. Tauris & Manch- 
ester (1998) examined this and found the fraction /(P) of 
pulsars beaming towards us in terms of pulsar spin period, 

/ (P) = 0.09 (log (P) - 1)^ + 0.03. (28) 

This equation shows that the faster the pulsar rotates the 
greater the chance that it is beamed towards us. The basic 
reason for this is that the pulsar cap angle (beam size or 
beam angle) is larger for fast rotators. A full explanation 
relies on details of the magnetic field and beam emission 
mechanism and we suggest the reader consider Tauris & 
Manchester (1998). However, the pulsar still may not be ob- 
served if selection effects conspire against our observing it. 
In fact there are some systems that we can rule out directly 
as not being observable in the radio regime. These are ac- 
creting systems or those in which the companions wind is 
energetic (lUarionov & Sunyaev 1975; Jahan-Miri & Bhat- 
tacharya 1994; Possenti et al. 1998). Thus, using Equation l28l 
we ignore those pulsar that are beamed away from us in 
combination with ignoring those that have a non-degenerate 
companion (including stars on the main sequence, although 
this may be too strict and could possibly be relaxed), or 
a pulsar which is accreting material. We also assume that 
it takes time for a pulsar to be observed (at radio wave- 
lengths) after the end of an accretion phase - here we as- 
sume it takes 1 Myr. Furthermore, any system which passed 
through a TZO phase is considered to be unobservable owing 
to the ejected nebulous material inhibiting radio transmis- 
sion. Model Fc results from these considerations. 

Because this model only culls the number of pulsars in 
the PP parameter space we do not show a comparison fig- 
ure. It is reassuring that the inclusion of beaming does not 
alter the good agreement found between the model and the 
observations in terms of the PP diagram. With the inclu- 
sion of beaming and the assumption that pulsars can not 
be observed for some time after accretion we find the bi- 
nary/isolated MSP ratio is 0.38 for Model Fc with the ma- 
jority of MSPs produced by wind accretion. 
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5.8.3 Propeller physics appraisal 

Here we examine further the effect propeUer physics (de- 
scribed in Section [3.3p has on the pulsar population. Moti- 
vation for this arises from the lack of binary MSPs coming 
from the RLOF channel in Model Fc. This arises because 
the propeller phase halts and then reverses the majority of 
spin-up phases (as seen in Figure In other words, many 
Roche-lobe overflow NS systems are having their spin-up 
phases truncated at some medium to slow spin period and 
therefore are not being 'observed' by binpop as they reside 
below the death line within the PP parameter space (at 
spin periods of ~ 0.03 or slower). This is in stark contrast 
to the standard method of assumed MSP formation, which is 
thought to occur via a low-mass X-ray binary phase (see Sec- 
tion 1 and references therein; although see also Deloye 2007 
for an interesting discussion on the low-mass X-ray binary- 
MSP connection) . This result begs the question 'what would 
these later models (that is, say. Model F) produce if pro- 
peller evolution was not included?'. To answer this we have 
repeated Model F with the propeller option turned off to 
create Model Fd. 

In comparison to Model F we find that the PP parame- 
ter space of Model Fd is similar, the most notable exception 
being a slight increase in the width of the pulsar bridge 
(in P). We certainly feel it is a plus to have the propeller 
option in our evolution algorithm, especially when one con- 
siders the claims of Winkler (2007), who suggests a new 
class of high-mass X-ray binary with the massive compan- 
ions losing mass in winds yet the NSs have spin periods of 
order 100—1000 s. Conversely, observations of X-ray binaries 
that contain a rapidly rotating NS, including the growing 
range of observed millisecond X-ray powered pulsars (Gal- 
loway 2007; Krimm et al. 2007), which we can produce in 
Model Fd, confirm that care must be taken when implement- 
ing propeller evolution. This suggests that improvements are 
required in how, and to what extent, propeller evolution is 
modelled. Our work presented here also corresponds to the 
claims made by Kulkarni & Romanova (2008). They show 
that unstable transfer of mass onto the accretor may occur 
- even during stable mass transfer events. Therefore, further 
investigation in this area is required, especially with regard 
to the possibility of transient pulsar mass accretion in X-ray 
binaries. 

It is interesting to note that with the increased number 
of RLOF binary MSPs in Model Fd that the binary/isolated 
ratio is 1.6. This compares well to the observed ratio (Huang 
& Becker 2007). However, it is too early on in our study to 
make conclusions based on such a comparison and there may 
well be a set of parameter changes that affect this result that 
we have yet to consider. 



5.8.4 Primordial parameter range 

Up until now the majority of our focus has been on pulsar 
evolution and how modifying the parameter values used in 
the assumed evolution of pulsars affects the distribution of 
pulsars that are produced. However, one of the advantages 
in coupling pulsar evolution with a full binary evolution al- 
gorithm is that it allows us to quantify the range of initial bi- 
nary parameters from which pulsars and particularly MSPs 
are generated. It also faciliates the investigate of how uncer- 



tain assumptions in binary evolution affect these ranges and 
the nature of the pulsar population. 

In Figure [T3] we show the distributions of initial binary 
parameters that lead to pulsar production in Models B, Fb 
and Fc. In these distributions we consider only those pul- 
sars that are plotted in their respective model PP diagrams 

- that is only those pulsars above the death line and for 
Model Fc those pulsars that satisfy the assumed selection 
effects. Also, a binary system is only considered once, even 
if that system produces two pulsars in the observable PP 
region (although this is relatively rare). We also include pul- 
sars that lie on the lower magnetic field limit. All three pri- 
mary mass distributions peak at around 8 Mq. This value 
is the minimum mass in which a star is guaranteed to form 
a NS (based on the stellar evolution algorithm). As we have 
addressed previously (Section 13. 5|) . the production of a NS 
from stars in the ~ 6 — 8 Mq primary mass range will depend 
upon the evolutionary pathway, particularly any associated 
mass loss/gain. The drop off at around 11 Mq in the mass 
distributions is chiefiy a stellar evolution feature - for solar 
metallicity models stars with M > 11 Mq ignite core-helium 
burning before the giant branch which gives a change in the 
radius evolution and stellar lifetime. Formation of NSs with 
primary masses less than 6 Mq can only occur via mass 
accretion or coalescence. The difference in the peak in and 
around 8 Mq between Model B and Fb is due to propeller 
evolution. Those NSs that accrete mass and go on to collapse 
further into BHs, in Model B, instead accrete less material 
due to the propeller phase in Model Fb and thus live longer 
lives as NSs and as such are more likely to be 'observed' 
within BINPOP. Although not examined in much detail here, 
at first glance one may consider this process to have a highly 
detrimental effect on the birth rate of BH low mass X-ray 
binaries when compared to NS low mass X-ray binaries (see 
Kiel & Hurley 2006 for discussion). However, because the NS 
is in a propeller phase for up to hundreds of millions of years 

- and will not be counted as a NS low mass X-ray binary 
during this time - these systems will not add greatly to the 
birth-rate calculations of these X-ray systems. In comparison 
to Models B and Fb, Model Fc producess a lower number of 
pulsars, as expected. The distribution of secondary masses 
for all three models peaks at around 5 Mq and covers the 
range from 0.1 — 33 Mq. This peak is due to a combination 
of the lower limit of 5 Mq for primary stars, the fact that 
dynamic or long lived mass transfer events occur most of- 
ten in binary systems in which the mass ratio is near unity 
and that the coelescence of two ~ 5 Mq stars leaves one 
that is able to form a NS. Finally we come to the orbital 
period distribution where we find that short orbital periods 
are favoured. Bearing in mind that the initial distribution of 
periods is fiat in log P, and thus, the relative number of bi- 
nary systems born with periods in the 1 — 10 d and 10 — 100 
d ranges are the same, it is interesting to see that the period 
distribution of pulsar progenitors does not follow that of all 
binaries. This indicates that close binary evolution is instru- 
mental in increasing the number of observable pulsars. With 
the advent of detailed pulsar evolution in combination with 
rapid binary evolution it is now possible (as demonstrated 
here) to consider the underlying main sequence stellar pop- 
ulation that generates the complete observed radio pulsar 
distribution. 
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Figure 13. Pulsar primordial parameter space distributions for 
primary mass, secondary mass and orbital period. The full squares 
represent the distribution of Model B, diamonds Model Fb and 
the plain line Model Fc. Distributions are normalised to the max- 
imum of each parameter across the three models. 



5. 8. 5 Observational constraints and predictions 

Conclusions on how complete our knowledge truly is for as- 
pects of binary evolution theory may be drawn if we com- 
pare certain pulsar population features directly with obser- 
vations. This may be completed in some limited form even 
though we do not yet consider selection effects fully. To do 
this we explore in greater detail the accretion physics at 
play in Model Fd. In particular, we examine MSP mass ac- 
cretion showing that under standard accretion physics as- 
sumptions the RLOF-induced MSPs accrete enough mate- 
rial for their magnetic fields to decay beyond observed values 
(in this model we assume a bottom field of 5 x lO'^ G). Fig- 
ure [14] shows that many MSPs that were spun-up via RLOF 
accreted greater than ~ 0.1 Mq worth of mass. Which, 
according to the prescription presented here is enough ac- 
creted material to decay/bury the field. In fact we find that 
~ 0.004 Mq is enough mass for a NS to accrete and sit on 
or close to the bottom field, or if we assume no bottom field 
it is enough mass that the NS will not be observable as a 
pulsar. We note that the MSP spin periods calculated in this 
work should be considered lower limits. The BSE algorithm 
takes care to conserve angular momentum and the spin pe- 
riod resulting from the accretion of material is based on the 
angular momentum transfered by the material. We also per- 
form a check that this does not violate energy conservation 
(which affects the lower left-hand region of Figure [H)) . how- 
ever, this check currently assumes that all the gravitational 
potential energy of the accreted mass is converted into rota- 
tional energy for the accreting star. Refining this treatment 
will be investigated further in the future. 

We also consider the mass range of pulsars produced in 
the above suggested model and compare this with observa- 
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Figure 14. Mass accreted since SNe for all pulsars rotating 
with spin periods faster than 0.02 s excluding those that fail our 
energy conservation check (see text for details). These pulsars 
are formed from Model Fd (with out propeller physics included). 
Black points represent those unobservable pulsars whose magnetic 
field has decayed to the assumed bottom magnetic field limit. All 
other points represent potentially observable MSPs. 

tions. Figure [15] depicts the distribution of NS mass against 
spin period and provides the number density of pulsars over 
the entire mass range. Note that in our models we assume a 
maximum NS mass of 3 Mq, which explains the upper cut- 
off in the mass distribution (and also affects how much mass 
a NS can accrete before becoming a BH). Overlaid on our 
model mass distribution is the observed pulsar masses with 
their errors as quoted from Table 2 of Nice (2006; and refer- 
ences therein) and Table 1 of Freire (2007; except for pulsar 
J1748-2021B). We also include the latest mass estimate of 
MSP J0437-4715 by Verbiest et al. (2007) and updates of 
J0705-hl807, J0621-hl002 and J1906+0746 from Nice, In- 
grid & Kasian (2008), while we make use of Lattimer & 
Prakash (2007) for pulsars J0045-7319 and J1811-1736. At 
this stage, due to the lack of observational constraints, we 
believe it would be inadvisable to put forward strong views 
on this matter. It is also instructive to see what occurs if 
we ignore those MSP-like NSs that would most likely be 
unobservable due to accretion induced field decay (again see 
Figure (15)) . We see that a large number of the highly recycled 
pulsars are now removed from the mass distribution. When 
only considering those potentially observable re-cycled NSs 
the apparent overlap of observed to theoretically produced 
pulsars, in the spin period ~ 0.003 to 0.01 range, is now 
more prominent. Beyond this spin range the number den- 
sity of observable model pulsars drops off somewhat. 

In keeping with our discussion on MSPs we now con- 
sider what MSP companion types are produced by Model Fd 
and the ratio of eccentric to circular orbits in MSP binaries. 
MSP binary systems comprise much of the interesting binary 
evolutionary phases in their histories and as such will play 
an important role in constraining binary and stellar evolu- 
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Figure 15. NS mass versus spin period for the complete popu- 
lation of Model Fd. Plus symbols represent NS formed from our 
model with Bs > 6 X 10^ G (and thus potentially observable) 
while the light dots represent all other model pulsars. The darker 
points show observations and their suggested errors. On the right 
we show the density of pulsar mass. We see that the majority of 
pulsars have masses in the range 1.3 to 1.5 Mq. 



tion theory via population synthesis. Figure \W\ depicts our 
Model Fd MSPs in the eccentricity-orbital period parame- 
ter space. The model MSPs in Figure[16]are selected to have 
Ss > 6 X 10^ G and spin periods less 0.02 s. Model Fd is 
also able to produce eccentric WD-MSP systems although 
the most eccentric of these systems in which the MSP is not 
in the process of accreting material is one with an eccentric- 
ity of 0.2. All eccentric MSP-MS binary systems depicted 
within Figure [16] are also accreting and thus not observable 
as radio pulsars. It is interesting to note the healthy amount 
of double NS systems Model Fd produces which include a 
millisecond pulsar and a large eccentricity. As an example to 
the extent of predictive power our future code will have on 
the possible understanding or interpretation of observations, 
we now consider the MSP J0514-4002A (Freire et al. 2004). 
This is the observed binary system depicted in Figure [TU] 
with an eccentricity of 0.888. At present this binary system 
has a measured mass of Msystcm ~ 2.453 Mq and resides 
within the globular cluster NGC 1851. The MSP, which ro- 
tates at P = 0.004991 s, has an estimated mass of < 1.5 AIq 
which leads Freire, Ransom & Gupta (2007) to argue that 
the minimum mass of the companion is 0.91 Mq. This be- 
ing the case Freire, Ransom & Gupta (2007) suggest that 
the companion is a heavy WD and that the most likely for- 
mation mechanism is via third body interaction due to the 
globular cluster high stellar density. Although the likelihood 
of this formation mechanism being the correct one is high 
they do suggest an alternative - a near head-on collision of 
a pulsar and giant star. However, because of where J0514- 
4002 A sits within Figure [TS] we consider a different possible 
binary companion type and formation mechanism for J0514- 
4002A. Some of the highly eccentric MSPs clustered around 
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Figure 16. MSP eccentricity versus orbital period for the MSP 
population of Model Fd. All pulsars here have magnetic fields such 
that Bs > 6 X 10'' G and spin periods less 0.02 s. Triangles repre- 
sent model double NS systems containing one MSP. Plus symbols 
are WD-MSP systems, symbols represent main sequence-MSP 
systems while solid sqares represents naked helium-MSP systems. 
There are not many naked helium-MSPs within this figure and 
all of them reside in circular orbits, thus rendering them hard to 
see. The black stars are observed pulsars taken from the ATNF 
pulsar catalogue (Manchester et al. 2005). 



1 < -Porb < 10 days within Figure [161 have low mass primary 
and secondary NSs, that is Msystcm < 2.7 Mq. The MSPs 
within this region have accreted only ~ 10~^ AIq and have 
been spun-up to periods, 2 x 10~^ < P < 10~^ s. So, we have 
some low mass systems which roughly fit the observed binary 
system containing MSP J0514-4002A, bearing in mind that 
our minimum initial NS mass is 1.3 Mq. This is especially 
the case if one assumes a minimum NS mass of ~ 1.2 Mq 
- the latest minimum pulsar mass estimate (Faulkner et al. 
2005). Prior to this it was assumed that pulsar masses were 
greater than 1.3 Mq - the value which BSE was calibrated 
to. This is something that will need to be revised in the fu- 
ture. Therefore, the model suggests that a low-mass highly 
eccentric MSP-NS system should also be considered as a 
possible explanation for this binary (in addition to the ec- 
centric MSP-WD scenario). We may also play this game for 
the Galactic disk eccentric millisecond binary system PSR 
J1903-I-0327 which has an orbital period of 95 days and an 
eccentricity of 0.44 (Champion et al. 2008). Although not 
plotted in Figure 1161 this system sits close to three dou- 
ble NS systems, perhaps suggesting that PSR J1903-I-0327 
could also be a low mass eccentric MSP-NS binary system. 
The model ratio of eccentric to circular MSP systems is 0.34. 
This is assuming that any system with an eccentricity less 
than 0.001 is counted as circular. The model does not follow 
eccentricity evolution below 0.001 owing to concerns about 
accuracy, i.e. systems with eccentricities lower than this are 
considered to be circular. 

Although not shown here directly it is possible for us to 



22 P.D. Kiel, J.R. Hurley, M. Bailes and J.R. Murray 



produce MSP-BH systems. Unfortunately, all of the MSPs 
in these systems accrete enough material that their mag- 
netic fields decay quickly beyond the pulsar visibility limit. 
However, it is interesting to see that we are able to produce 
rapidly rotating NSs with BH companions. A typical evolu- 
tionary formation scenario begins with two large stars, the 
primary Mi ~ 65 Mq and secondary M2 ~ 30 Mq in a close 
orbit separation a ~ 40 -Rq (an orbital period of ~ 5 days). 
Both stars begin to rapidly lose mass in winds and eventu- 
ally the primary inititates mass transfer, at a system time 
of 2.5 Myr by which time it has already lost ~ 4 Mq and 
the companion has lost ~ 0.3 Mq. When the RLOF mass- 
transfer phase ends some 2 Myr later the primary star is now 
a core-helium burning star with 8 Mq and the secondary is a 
66 Mq main sequence star. The orbit has increased to a sep- 
aration of ~ 180 -Rq. The primary star then rapidly evolves 
through naked helium phases to become a 2.4 Mq NS. The 
system survives the SN event and is now 5.2 Myr old with a 
separation of ~ 400 -Rq between the stars. The companion 
begins to evolve quickly off the main sequence and in the 
0.5 Myrs it takes to become a 13 Mq BH it loses 30 Mq 
in a wind. The NS has accreted 0.07 Mq of this material 
which is enough to spin it up to a millisecond spin period. 
We find that 0.8% of the binary pulsars rotating with pe- 
riods less than 0.02 have a BH companion (including those 
that have decayed to the bottom magnetic field limit). Of 
these rapidly rotating NSs with BH companions just under 
half occur with orbital periods less than 50 days (the small- 
est orbital period is just under 0.1 days). Also, the majority 
of them (94%) have eccentricities greater than 0.2. It can 
be supposed that if one modifies our method for decaying 
the magnetic field during accretion such that (some of) the 
above system MSPs do not decay beyond the pulsar death 
line, it could be possible to produce observable MSP-BH 
systems in relativisitic orbits (that is a heavy, tight and ec- 
centric binary system in which one may constrain general 
relativity) within the Galactic disk. 



6 DISCUSSION 

In the previous sections we investigated a range of input pa- 
rameters and evolutionary assumptions that affect the out- 
comes for isolated and binary pulsar evolution. This lead us 
to a preferred model which gave good coverage of the PP pa- 
rameter space. For that model we examined NS and orbital 
properties comparing them to observations. We next inspect 
some aspects of the results in further detail and comment 
on areas we think are important for future work. 

Modelling beaming and accretion selection effects signif- 
icantly reduces the number of observable pulsars predicted 
by our model. Within this Model Fc there was uncertainty in 
how to deal with Thorne-Zytkow objects. Within Section[5]6] 
our assumption was that the NS reset during the destruction 
of the TZO (made for simplicity) . At present we are unsure 
if this is a realistic assumption to make. If we instead assume 
the pulsar is unchanged when passing through a TZO evo- 
lutionary phase we find many more isolated MSP systems 
are produced. Coupled to this, in Section [5]8]2] we assumed 
that when the material is expelled from the TZO, leaving be- 
hind a NS, that the expelled material enshrouds the NS for 
long enough that it is not observed before crossing the death 



line (very much like the assumption for post-accretion NSs) . 
Again, we are unsure how physically correct this assump- 
tion is. Below are a number of other possible evolutionary 
assumptions that could be made when dealing with these 
systems. Perhaps the NS should accrete some amount of ma- 
terial, decaying the magnetic field, and, if enough material 
is accreted, form a BH (Fryer, Benz & Herant 1996). Or if a 
BH is not formed perhaps a revitalised NS is born (Podsiad- 
lowski 1996) - a standard isolated pulsar, or a pulsar which is 
below the death line. But one may even consider an outcome 
that is quite interesting. That is the possible formation of 
a magnetar. It could even be assumed that the surrounding 
material causes a severe braking mechanism, the outcome 
of which is a young hot highly magnetic slowly-rotating NS 
(see below for further discussion on possible magnetar pop- 
ulation synthesis). 

We showed in Section 15.51 that values of the magnetic 
field decay timescale, re, below 100 Myrs do not provide a 
realistic distribution of pulsars. Then in Section [5.61 we ex- 
plained that using tb = 1000 Myr does not cover enough 
of the slow rotating pulsars in the PP diagram but that 
using Tb > 2000 Myr over corrected for this. For our pre- 
ferred model we found tb = 2000 Myr to be optimal. We 
note that the choice of maximum age for the population, in 
our case 10 Gyr, affects the time available for decay and as 
a result affects the appearance of the pulsar distribution at 
the lower edge of the main pulsar island. However, the ap- 
pearance of this region also depends crucially on the details 
and/or shortcomings of the assumed death line used in the 
model (which we now discuss). 

In all our models the recycled pulsars are not suitably 
limited by the death line. Of course, the inverse Compton 
scattered death line models have not been considered in any 
detail, and as Hibschman & Arons (2001) suggest all three 
models may be required to suitably simulate the complete 
pulsar population death line. In fact, Hibschman & Arons 
(2001) show regions of PP in which the different methods 
would come into play. The standard pulsar island is shown 
to be modelled by CR, the recycled pulsars by nonresonant 
inverse Compton scatter and the high magnetic field pul- 
sars governed by the resonant inverse Compton scattered 
method. This matches our findings that CR models the stan- 
dard pulsar island death line in a suitable manner but is too 
lax in culling recycled pulsars. To test this further we now 
examine what effect decreasing /^'m to 0.10 (from the pre- 
viously assumed 0.15) has on our pulsar population. The 
recycled pulsar island is now more effectively constrained 
by the death line. However, the death line is now too effec- 
tive at culling pulsars from the standard pulsar island. As 
explained previously (see section |377|| . the death line is bro- 
ken up into two regimes, saturated and unsaturated. Within 
all the models so far it seems that the majority of pulsars 
are constrained by the saturated death line region (even if 
we are to decrease /^IIJi to 0.03). In fact, any pulsar at the 
begining of its accretion phase is governed by the saturated 
regime (although at some point after this it may move into 
the unsaturated regime for a time). To increase the infiuence 
of the unsaturated regime we can take some liberties with 
the assumed unsaturated /^Im value and let it approach 
zero. However, this eliminates a wedge of pulsars from a re- 
gion of the PP space in which pulsars are observed and this 
is not an acceptable parameter change. We note that death 
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line calculations are model dependent and in particular the 
assumed NS equation of state is an important feature and is 
something we have not considered here in any detail. Further 
work is required in examining all the Harding, Muslimov & 
Zhang (2002) death lines but is beyond the scope of this 
paper. 

All the models presented here allow the creation of NSs 
with spin frequencies in excess of 1000 Hz. The only reason 
these NSs do not spin beyond ~ 3300Hz (see, say. Figure [10)) 
is that, within BSE, no star is allowed to spin beyond their 
break-up speed. There are theoretical mechanisms that al- 
low the regulation of pulsar spin and restrict it to some 
spin velocity less than the break-up velocity. One such pos- 
sible mechanism is gravitational emission from the pulsar 
itself. There are two theories as to how a pulsar may emit 
gravitational radiation and thus self-regulate its own spin 
period. One theory suggests that if a pulsar surface is not 
smooth or purely spherical (this would be the case for those 
that have accreted material or that have strong magnetic 
fields and are rotating rapidly) then this Eisymmetry in the 
pulsars features may be mapped to the pulsars space time 
curvature and cause energy disipation in the form of gravity 
waves (Shapiro & Teukolsky 1983; Bildsten 1998; Payne & 
Melatos 2005). The other theory assumes that the Rosby 
wave (r-wave or r-mode) instability (Levin 1999; Friedman 
& Lockitch 2001) is not damped. This instability basically 
arises due to the form of the wave and the speed of rotation 
of the pulsar. Chandrasekhar (1970) showed that Maclaurin 
spheriod^ with uniform density, rotating uniformly are un- 
stable to 1 = m = 2 polar modes (i.e., those waves in which 
the wave front stretches from pole to pole and travels around 
the sphere). In a non-rotating or slow rotating star, a par- 
ticular mode, the non-axisymmetric instability removes pos- 
itive and negative angular momentum from the forward and 
backward mode respectively (Friedman & Lockitch, 2001), 
thus all the non-axisymetric modes cancel out and no insta- 
bility forms. In a rapidly rotating star (~ that around the 
maximum observed pulsar spin rate; Levin 1999), a back- 
ward moving mode may be moving forward as observed by 
an inertial observer and therefore radiate away positive an- 
gular momentum - in effect increasing the mode amplitude 
(Friedman & Lockitch 2001). In this manner gravitational 
radiation drives the mode and the star spins down until the 
instability dies (Levin 1999). 

With the suggested discovery of a pulsar with sub- 
millisecond period (Kaaret et al. 2007), the maximum spin 
rate possible for a pulsar is work we would like to consider in 
the future. At the forefront of this work will be contempla- 
tion of the above gravitational radiation pulsar spin-down 
methods. However, we wish to point out that the assumed 
accretion-induced magnetic field decay parameter may play 
an important role in determining the observed MSP spin 
period cut off (see Section 15. 3p and is something that will 
need to be considered in future work. 

In recent years magnetars have become a distinct sub- 
class of the pulsar population. These systems are highly 



^ A method for modeling uniform density ellipsoids. In particular 
it gives a description of the relationship between the eccentricity 
of the system, in this case the pulsar, and its uniform angular 
momentum. 



magnetic NSs with magnetic fields of order 10 — 10 G. 
Magnetar formation is not well understood, though work 
by Thompson & Duncan (1993) and Duncan & Thompson 
(1992) showed that nascent NSs which rotate rapidly (< 0.03 
s) can have magnetic fields within the above range. In terms 
of the work being completed here we have not yet consid- 
ered the possibility of modeling these highly magnetised, 
slow rotating pulsars. The reasoning behind this is that 
we are trying to produce radio pulsars and these systems 
have only been found in X-rays (anomalous X-ray repeaters) 
and gamma-rays (soft gamma-ray repeaters: see Woods & 
Thompson 2006 for review). However, in light of the now 
constant flow of new and improving theoretical and obser- 
vational considerations of these objects (Spruit & Phinney 
1998; Heger et al. 2003; Heger, Woosley & Spruit 2005; Mez- 
zacappa 2005; Thompson 2006; Harding & Lai 2006; Blondin 
& Mezzacappa 2006) an attempt to model these systems 
with very simple assumptions is now possible and will be 
included in future work. 



7 CONCLUSIONS 

We model the evolution of radio pulsars, both single and 
in binary systems. Models are evolved from realistic initial 
mass functions, period distribution and parameter ranges. 
Our primary goal is to demonstrate that we have an evo- 
lution algorithm that can generate the range of pulsars ob- 
served in the Galaxy, thus laying the groundwork for a syn- 
thetic Galactic population code that will ultimately include 
kinematics and selection effects. However, through inves- 
tigation of the PP diagram (comparing model-model and 
model-observation) we have already been able to demon- 
strate how uncertain factors of pulsar and binary evolution 
can be constrained. We have also been able to make prelim- 
inary predictions about the Galactic pulsar population, in 
particular MSPs. This represent the first study combining 
detailed binary evolution and pulsar physics. Of course a 
study such as this is hampered by many poorly constrained 
parameters which may result in misleading results. This is 
especially true if the complete relevant parameter space is 
not examined and multiple observations are not used as 
benchmarks. Therefore, we once again remind the reader 
that some caution is required when examining our results. 

In terms of pulsar evolution it has been possible to pro- 
vide constraints on various aspects of the modelling. For in- 
stance, our models can demonstrate that magnetic field de- 
cay timescales of order 100 Myr or less can not produce the 
required distribution of spun-up pulsars. We also show that 
propeller evolution has a large impact on pulsar populations 
and can be an important component of pulsar modelling. In 
particular we show that the inclusion of the propeller mech- 
anism causes wind accretion to be the dominant pathway 
for MSP production. In fact, neglecting propeller evolution 
is necessary in order to produce any sort of RLOF disk ac- 
cretion powered MSP (the typical type of observed accre- 
tion powered MSP). To further investigate this aspect there 
is a need to examine methods for dealing with accretion 
physics in more detail - including the addition of transient 
RLOF accretion features. We find that accretion-induced 
magnetic field decay (or apparent field decay) provides a 
natural method for forming rapidly rotating low magnetic 
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field NSs. Without any such magnetic field decay we find 
there is no hope of theoretically reproducing the observed 
pulsar PP distribution. The accretion-induced method as- 
sumed here works well for all models performed and we de- 
pict how changes in the officcncy of the decay modifies the 
populated PP region. Including an angular momentum wind 
accretion efflcency parameter - lowering the angular mo- 
mentum accreted - assists in producing some of the tightest 
observational PP constraints. We also find that the addi- 
tion of an electron capture assymctric SN kick distribution 
results in the formation of a greater number of MSPs than in 
previous models. As such this evolutionary assumption will 
doubtlessly be an important feature in future works and re- 
quires further examination. 

We demonstrate the predictive power of our code by 
examining the pulsar mass range and eccentricity of MSP 
systems produced in a favoured model. We must wait to im- 
pose selection effects before we can advise as to what frac- 
tion of these arc potentially observable. We also show that 
under standard accretion physics assumptions NSs are able 
to accrete up to ~ 1 M© of material, which can spin them 
up to rapid rotation rates and decay their magnetic fields to 
strengths considered insufficient for producing radio emis- 
sion. 

We conduct an initial investigation into how the PP 
parameter space is modified when selection effects are con- 
sidered. These include beaming and blanketing of the pulse 
by material thrown off by stellar winds or during mass ac- 
cretion events. We study the distributions of primordial pri- 
mary mass, secondaxy mass and orbital period for binaries 
that produce pulsars. Changes in these distributions due to 
assumed evolutionary variations are examined. In particular 
we show that close binary evolution is effective in boosting 
the numbers of observable pulsars. Finally, it is clear from 
all of our models that a greater level of detail is required 
when calculating the complete pulsar population death line. 

The work shown here makes use of the first of three 
modules that will eventually be able to simulate pulsar ob- 
servational surveys, including next generation radio surveys 
such as the Square Kilometer Array. This involves following 
the positions of pulsars within the Galactic gravitational po- 
tential and modelling survey selection effects. The long term 
aim is to include multi-wavelength survey capability into our 
models, and thus observe any Galactic stellar population. It 
is important to point out here that once complete this code 
will consist of the most comprehensive treatment of input 
physics and selection effects for pulsar simulation, something 
not previously attempted. 
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